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ABSTRACT 


Two  variants  of  the  multivalue  method,  the  correc- 

M 

tor  only  and  the  P  (EC)  methods,  for  the  numerical 
solution  of  the  differential  equation 

y(p)  =  f (t,y , . . . ,y (p“q) ) 

are  studied.  It  is  shown  that  the  numerical  solution  of 
a  convergent  corrector  only  method  is  also  the  solution 
of  a  stable  multistep  method.  A  useful  technique  for 
analyzing  multivalue  methods  is  discovered  which  makes 
it  possible  to  formulate  a  constructive  definition  for 
the  order  of  a  method.  When  the  stepsize  is  fixed, 
necessary  and  sufficient  conditions  for  rth  order  q- 
convergence  are  given,  and  explicit  expressions  for  the 
asymptotic  behaviour  of  the  global  error  are  derived. 
When  the  stepsize  is  varied  in  a  reasonable  way,  it  is 
shown  that  the  property  of  rth  order  q-convergence  is 
preserved  for  q  =  1  but  not  for  q  >  1.  In  the  special 
case  of  a(k+l)-value  Adams  method,  stability  can  be 
guaranteed  by  keeping  the  stepsize  constant  for  at  least 
k-1  steps  after  a  change  in  stepsize. 
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CHAPTER  1 


INTRODUCTION 

The  objective  of  this  thesis  is  to  study  multi¬ 
value  methods  for  solving  ordinary  differential  equa¬ 
tions  of  the  form 

y(p)  =  f(t,y . . (P-<3) >  . 

Multivalue  methods  have  been  developed  by  Gear  (1971) 
as  generalizations  of  Nordsieck' s  method,  which  saves 
approximations  of  the  higher  order  derivatives  of  the 
solution  without  directly  evaluating  these  derivatives. 
In  his  original  paper  Nordsieck  (1962)  demonstrates  the 
equivalence  of  his  method  to  an  Adams-Moulton  method. 
Subsequently  Osborne  (1966)  showed  that  the  Nordsieck 
methods  were  equivalent  to  multistep  methods. 

Gear  (1967)  ,  in  extending  the  work  of  Descloux, 
noted  that  any  multistep  method  can  be  written  in  many 
forms.  To  take  a  step  requires  a  set  of  back  values 
and  derivatives,  and  this  set  defines  an  interpolating 
polynomial.  There  are  many  ways  to  represent  this 
polynomial,  and  in  particular,  Nordsieck’s  method  saves 
the  higher  order  derivatives  of  the  polynomial  evaluated 
at  the  latest  value  of  the  independent  variable  t. 
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These  saved  values  are  denoted  by  the  vector 


-n  =  C0l(yn'hyA . ^  lynk+P  1,/(k+p-l)!) 


where  h  is  the  stepsize.  This  particular  formulation 
of  a  (k+p) -value  method  is  useful  because  it  makes 
stepsize  changing  easy. 

Gear  (1971)  has  given  the  theory  of  multivalue 
methods  considerable  attention  but  some  of  the  analysis 
remains  incomplete.  In  particular,  the  classical  result 
of  "stability  and  consistency  are  equivalent  to  conver¬ 
gence"  has  not  been  obtained. 

The  major  problem  in  obtaining  this  theorem  is 
that  of  formulating  an  adequate  definition  of  consis¬ 
tency.  In  Gear's  book  there  are  two  definitions  for 
rth  order  consistency,  one  for  the  case  r  <  k  and  another 
for  the  case  r  >  k.  For  r  >  k  the  definition  of  rth 
order  consistency  is  difficult  to  apply  because  it  is 
not  constructive.  And  for  r  <  k  the  definition  of  rth 
order  consistency  is  not  satisfied  by  all  methods  having 
convergence  of  order  r.  As  an  example  consider  the 
following  multivalue  method  for  first  order  equations: 

.  ,2 

hi  .  h  » 

yn  “  yn-l  +  2  yn-l  +  T  yn-l  ' 


hyn=hf(tn'yn> 


t 


■ 


3 


hf  (t 


n'V 


The  zeroth  component  of  the  local  truncation  error  x 

— n 

is  defined  by 


y(V  +  (V0  =y<tn-l)  +  !  y  (tn-l)  +  T-  y"(tn-l)' 


and  since  (t  )  is  not  of  order  h  .  the  method  is  not 
-n  o 

consistent  in  the  usual  sense.  With  a  change  of  nota¬ 
tion  the  method  can  be  written 


yn-l 


I 


hy  ’  =  hf (t 

n 


/ 


hy ' =  hf (t 
Jn 


n'V 


This  example  makes  it  clear  that  the  method  is  equiva¬ 
lent  to  Euler's  method  with  a  starting  error  of  order 
h,  and  therefore  the  method  is  convergent. 

In  this  chapter  q-convergence  is  defined  for 
multivalue  methods.  The  formula  used  by  a  multivalue 
method  is  written  as  the  sum  of  a  linear  part  plus  a 
nonlinear  part.  When  f  =  f  (t) ,  the  formula  becomes  a 
linear  difference  equation  having  a  unique  solution  of 
a  special  form.  The  order  of  the  linear  part  of  the 
formula  is  defined  to  be  the  extent  to  which  the  solution 


of  the  difference  equation  agrees  with  the  correct 

solution.  Also  in  this  chapter  it  is  shown  that  the 

numerical  solution  of  a  multivalue  method  is  the 

solution  of  a  corresponding  multistep  method. 

In  Chapter  2  the  q-root  condition  and  the  rth 

q-order  condition  are  defined,  and  it  is  proved  for 

M 

both  the  corrector  only  and  the  P  (EC)  methods  that 


r — 1 

q-convergence  like  o (h  ) 


q-root  condition 

rth  q-order  condition  . 


In  Chapter  3  it  is  shown  that 


q-root  condition 


* 


q-convergence  like  0 (h  ) . 


rth  q-order  condition 


Also  in  Chapter  2  an  example  is  given  of  a  weakly  stable 
(p+2) -value  method  for  pth  order  equations  which  is  1- 
convergent  of  order  4. 

In  Chapter  4  the  asymptotic  form  of  the  error  as 
h  0  is  obtained  for  a  special  class  of  starting  values. 
This  result  is  applied  to  the  case  of  correct  starting 
values  for  strongly  stable  methods.  Partial  results  in 
this  area  are  obtained  by  Gear. 
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Often  in  practice  the  stepsize  is  varied.  Let 

hn  =  Then  the  interpolatory  technique  for 

varying  stepsize  is  accomplished  by  premultiplying  ^ 

by  the  matrix  C  =  diag(l,h  /h  (h  /h  )k+P_1) 

n  3  n'  n-1'  n'  n-1 

before  taking  the  next  step.  If  a  method  is  1-conver- 
gent  of  order  r  for  fixed  stepsize  then  it  is  also  1- 
convergent  when  the  stepsize  changes  satisfy  the  follow¬ 
ing  restrictions: 

1.  h  >  Ah 

n 

N 

2.  y  In  -h  , I  <  Vh 

L ,  1  n  n-11  - 

n=l 

where  h  is  the  maximum  stepsize,  A  and  V  are  fixed  posi¬ 
tive  constants,  and  N  is  the  number  of  steps.  For  q>  1, 
q-convergence  of  order  r  is  generally  not  preserved 
when  the  stepsize  is  varied.  In  fact,  almost  any  method 
which  satisfies  the  definition  of  rth  order  q-convergence 
for  variable  stepsize  is  also  1-convergent  of  order  r. 

In  the  second  section  of  Chapter  5,  it  is  shown  that  a 
(k+1) -value  Adams-Bashforth-Moulton  method  for  first 
order  equations  is  stable  if  the  stepsize  does  not 
change  more  often  than  every  k-1  steps.  The  results 
of  Chapter  5  are  improvements  of  work  done  by  Tu  (1972) . 


;3taj 
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Section  1.1  The  Problem 

We  are  interested  in  constructing  an  approximate 
solution  for  an  initial  value  problem  which  conforms  to 
the 


Definition  of  Lipschitz  problem: 

Let  p  and  q  be  fixed  positive  integers  where 
q  <  p.  The  pth  order  q-dif f erential  equation  y^  = 
f (t,y , . . . ,y )  is  called  a  Lipschitz  differential 
equation  if 

1.  f  is  a  continuous  function  of  t  on  [0,1]. 

2.  f  is  uniformly  Lipschitz  continuous  in  the 

/  •  \ 

y  's;  that  is,  there  exist  constants  Lq, 

L,  , . .  .  ,L  such  that 

1  p-q 


~  (i) 

f  (t  ,y ,  .  .  .  ,y  , 


fy(p  q) ) _  f  (t,y,...,y(l)  ,...,y(p~q))  | 


„  T  |  ~ (i)  (i) 

<  L^ | y  -  y 


for  all  y , . 


~(i)  (i)  v(p-q) 

•  tY  'Y  r  •  •  •  rY 


and  for  all  .0  <  t  <  1. 


A  Lipschitz  problem  consists  of  a  Lipschitz  differential 
equation  together  with  p  initial  values  Yq  'Yq  /  •  •  •  'Yq^  ^  .  □ 


For  notational  simplicity,  we  consider  only  one 
differential  equation.  The  theory  extends  almost  imme¬ 
diately  to  systems  of  equations  if  absolute  values  are 


■ 


. 


are 
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replaced  by  norms  and  derivatives  like  8f/Sy^ 
considered  to  be  Jacobian  matrices. 


Section  1.2  Polynomial  Methods 

Polynomial  methods  are  discussed  in  order  to  mo¬ 
tivate  the  definition  of  multivalue  methods  and  the 
definition  of  q-convergence . 

Let  the  solution  to  the  Lipschitz  problem  be 
y (t) .  Our  aim  is  to  construct  approximations  y1(t;h) 
to  y  ^  (t)  on  the  interval  [0,1]  for  i  =  0,l,...,p-q. 
Higher  derivatives  of  y  are  not  approximated  because 
they  are  not  required  for  the  evaluation  of  f. 

First  let  us  introduce  the  grid  of  equidistant 
points  0  "  tg  <  t^  <  . . .  <  tN  =  1  and  define  the  stepsize 
h  =  1/N.  The  approximate  solution  y1(t;h)  is  constructed 
by  piecing  together  N+l  polynomials  each  of  degree 
<  k+p-1  where  k  >  1: 

yi(t;h)=  p^^  (t)  for  tn~  j  <  t  <  tn+^  . 

Three  variants  of  the  polynomial  method  are  defined  which 

M 

are  analogues  of  multistep  corrector  only,  P{EC)  ,  and 
M 

P (EC)  E  methods.  All  three  methods  successively  deter¬ 
mine  a  sequence  of  polynomials  p^  (t) ,  p^  (t ) , . . . ,pN  (t) . 

But  first  Pg  (t)  must  be  obtained  by  some  starting  pro¬ 


cedure  . 
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In  order  to  justify  our  definition  of  convergence, 
let  us  write  the  differential  equation  as 


(P  D  f4-\-  f(T,y(T), - ,y  (T)  )dx  . 

/ 

0 


y  '*  '  (t)  = 


Then  it  is  not  necessary  for  f  to  be  a  continuous  function 
of  its  first  argument.  At  worst,  y^  ^  (t)  might  only 
be  absolutely  continuous  in  which  case  even  for  the  best 
polynomial  approximation. 


max 

0  <  t  £  h/2 


P0  (t)-y (t) 


o(hp_1) 


is  the  strongest  statement  that  can  be  made. 


Definition  of  q-convergence : 

A  polynomial  method  is  q-convergent  if  for  any 
Lipschitz  problem, 

y^  (t ;h)  =  y^  (t)  uniformly  in  t 

for  i  =  0,1,..., p-q  whenever  pQ  (t)  is  a  function  of  h 
satisfying 

P0  (t)  =  y(t)  +  o(hp_1)  uniformly  on  [0,h/2].  □ 


Definition  of  rth  order  q-convergence  for  polynomial 

methods : 

Let  r  be  a  positive  integer.  A  polynomial  method 
is  q-convergent  of  order  r  if  for  any  Lipschitz  problem 
with  solution  y  (t) eCr+^ [ 0 , 1] , 

y^(t;h)  =  y^  (t)  +  0  (hr)  uniformly  in  t 

for  i  =  0,1,..., p-q  whenever  Pq (t)  is  a  function  of  h 
satisfying 

Pg  (t)  =  y(t)+0(hr+p  uniformly  on  [0,h/2].  □ 

Throughout  this  thesis  the  use  of  the  little-oh  or  the 

big-oh  notation  carries  with  it  the  assumption  that  the 

limit  or  the  bound  is  uniform  with  respect  to  time 

(t  and  n) .  Nevertheless  there  may  be  a  dependence  on 

the  method  and  the  problem. 

The  accuracy  of  the  approximate  solution  is 

limited  by  the  use  of  polynomials  of  degree  k+p-1. 

The  initial  value  y^  ^  (t)  is  approximated  on  [0,h/2] 

by  a  polynomial  of  degree  k,  and  the  error  in  such  an 

k+1 

approximation  is  normally  of  order  h  .  Therefore  no 
polynomial  method  has  order  of  q-convergence  greater 


than  k+1 . 
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The  corrector  only  variant  of  the  polynomial 
method  defines  successive  polynomials  by 

t-t 

p_  (t)  =  p  ,  (t)  +  to  L 
n  n-l  n  h 

where 

L(s)  E  £k+p-l  sk  P  1  +  • • •  +  £is  +  £0 

and  to  is  chosen  so  that  p  (t)  satisfies  the  differen- 
n  n 

tial  equation  at  t  =  t  . 

Consider  the  equation  y^  =  f  (t)  ,  and  suppose 
that  p  ^(t)  is  known.  To  determine  pn(t),  is  chosen 

so  that  the  differential  equation  is  satisfied  at  t  =  t  : 


f  (t  )  =  p^^  (t  )  =  p^l  (t  )  +  to  Z 

v  n  ^n  v  n  *n-lv  n  n  up  p 


In  general,  Pj[?}(tn)  ?  f  (tn)  ,  and  so  must  not  be  zero. 

Since  multiplying  L(s)  by  a  nonzero  constant  does  not 

change  the  method,  we  assume  that  L(s)  is  normalized  so 

that  Z  =  1 . 

P 

The  polynomial  L(s)  determines  what  information 
is  retained  when  one  step  is  taken.  To  be  more  specific, 
the  ith  derivative  of  p  1(t)  is  equal  to  the  ith  deriva¬ 
tive  of  p  (t)  at  the  roots  of  ((t-t  )/h). 

n  n 
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As  an  example,  consider  the  explicit  k-step 
Adams-Bashforth  method,  which  uses  hy '  ...  hy ' 

^n-k '  and  ^n-1  ^e^ermine  Yn  •  The  saved  values 
define  an  interpolating  polynomial  Pn-1  (t)  so  that 

pn-l(tn-l>  =  yn-l  and  Cl  (tn-i>  =  C-i  for  1  =  1 '2 . k- 

For  explicit  methods,  y^  is  obtained  by  extrapolation; 
that  is,  yn  =  P„_1(tn).  Hence  pn(tn>  =  p„_1 (tn) ,  and 
so  L(0)  =  0.  Also,  since  pn(t)  is  determined  by  hy^, 

hyn-l'  *  *  *  'hyn-k+l'  and  Yn'  it:  follows  that  Pn(t)=  P^-i 
for  t  —  t^_^  r  2  t  •  •  •  , Pj^—k+l *  Thus  L'(s)  —  0  for  s  = 

“1 ,  -2 ,...,- (k  -  1) .  These  k  conditions  are  sufficient 
to  determine  a  normalized  polynomial  L(s)  of  degree  k. 

As  a  second  example,  the  implicit  k-step  Adams - 

Moulton  method  uses  hy'  ,,  hy'  0,...,hy'  .  ,  and  y  n 

■“n-r  Jn-2'  '  ^n-k'  in-l 

to  determine  y  .  Just  as  in  the  first  example,  L' (s)  =0 
for  s  =  -1,  -2 (k-1)  .  Only  one  more  condition  is 
needed  to  determine  L(s),  and  that  is  to  choose  L(0)  so 
that  the  formula  is  of  order  k+1. 

One  way  of  representing  the  polynomial  pR  (t)  has 
been  illustrated.  Another  more  convenient  representa¬ 
tion  for  p  (t)  is  the  vector  of  scaled  derivatives 
n 

,  k+p-1  ,,  ,  » 

=  ool(pn(tn,'hpA(tn) .  (k+p-1)!  p„  P’  (tn})  ‘ 


The  method  is  said  to  be  in  normal  form  if  the  polynomial 
is  represented  by  scaled  derivatives. 


2) 


Let  us  write  the  method  in  normal  form: 


(1.1)  a  =  Aa  n  +  03  l 

— n  — n-1  n— 

where  A  is  the  (k+p)x(k+p)  Pascal  triangle  matrix 


'  1  1  1  . 
1  2  . 

1  . 


.  k+p-1 

'k+p-r 

e 

l  2  J 


\ 


0 


1 


and  e  col  (£Q  ,  r  .  .  .  ,  .  (Indexing  starts  from 

0  for  all  vectors  and  matrices.)  Since  03n  is  chosen  so 
that  Pn(t)  satisfies  the  differential  equation  at  t  =  t  , 


2-^  (a  )  =  f(t  ,(a  )  ,  .  .  .  ,  ^  (a)  )  . 


hP  ~n  p 


n'  — n 


0 


p-q  -n  p-q 


Therefore 


0)  = 
n 


hP 

£7-  F  ( t,  A  a  n  +  o)£) 
p!  n  — n-1  n— 


(Aa  ,  ) 
—n-1  p 


where 


F ( t ,  a )  =  f  (t,aQ , .  .  .  , 


(p-q) 1 


a  )  . 


D-q  p-q 
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A  unique  solution  exists  for  03n  if  h  is  sufficiently 
small  so  that 


(1.2) 


p-q 

I 

i=0 


UjL.hP 


“I 


<  1 


Note  that  if  ln  =  =  .  .  .  =  Z  =  0  then  co  is  deter- 

0  1  p-q  n 

mined  explicitly.  Thus  there  are  explicit  and  implicit 
methods . 

For  implicit  corrector  only  methods ,  o>  is  often 
the  solution  of  a  nonlinear  equation,  which  must  be 
solved  iteratively.  Predictor-corrector  polynomial 
methods  are  developed  for  this  purpose.  For  values  of 
h  satisfying  (1.2),  the  following  iteration  converges: 


0) 


n,  (0) 


=  0  , 


oj  /  ,  n  x  =  — r  F  (t  ,A  a  ,  +  03  ,  x  i)  -  (Aa  n  ) 

n,  (m+1)  p!  n'  -n-1  n,(m)—  — n-l  p 


Define  a  ,  x  =  Aa  ,+03  ,  ,  and  then  the  corrector 

— n,  (m)  — n-1  n,(m)  — 

iteration  becomes 


d.3a)  an,(0)  n-1 


d.3b)  — n ,  (m+1  )  — n,  (m)  +  —  [p!  F(tn'^n,(m))  (-n,(m))p] 


In  practice  only  a  finite  number  of  corrector  iterations 


■ 
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are  done,  and  to  simplify  the  analysis,  assume  that  the 
number  of  iterations  is  a  constant  M: 


(1.3c) 


a 

— n 


^n,  (M) 


These  last  three  equations  constitute  the  definition  of 
a  P  (EC)  method.  The  abbreviation  stands  for  "predic¬ 
tion,  (evaluation  of  f,  correction)^. " 

Note  that  in  general,  the  polynomial  represented 
by  does  not  satisfy  the  differential  equation  at 

t  =  t  .  This  can  be  remedied  by  performing  a  final 
correction  on  the  pth  component  of  ^ : 

hp 

^n  “  ^n , (M ) +  — p  F (tn'-n, (M) } ~ (-n, (M) 


where  6_  is  a  column  vector  with  a  1  in  the  pth  position 
P 

(zero-origin  indexing)  and  zeros  elsewhere.  This  variant 

M 

of  the  polynomial  method,  called  a  P  (EC)  E  method,  is 
not  investigated  because  of  the  extra  complications 

M 

involved.  The  analysis  is  similar  to  that  for  P (EC) 
methods . 


Section  1.3  Multivalue  Methods 

Multivalue  methods  are  generalizations  of  poly¬ 
nomial  methods  in  normal  form.  A  multivalue  method  is 
defined  either  by  (1.1)  or  by  (1.3)  where  A  is  replaced 


' 
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by  the  general  matrix  A  and  for  a  P(EC)M  method.  £  is 

P 

not  required  to  be  1.  Methods  with  £  ^  1  are  less 

P 

accurate  and  less  stable,  and  thus  they  are  less  likely 

to  have  practical  value.  Both  Gear  (1971)  and  Stetter 

(1973)  consider  the  possibility  that  Z  ?  1 ,  and  so  for 

P 

completeness,  we  do  likewise. 

The  numerical  solution  produced  by  a  corrector 
M 

only  and  a  P  (EC)  method  satisfy  respectively 


a  =  (I  -  Z  6T) Aa  ,  .  „ 

— n  - p  — n-1  —  p 


hp 

+  Z  F  (t 


,  a 
n  — n 


and 


-n 


T  M  t  M-m 

=  (I-  £  61)  Aa  ,  +  l  (I-  Z  6„)  Z  x 

- p  -n-1  m^1  - p  -  pi 


x  F 


(t  ,  a 
n  '-n 


(m-1) 


) 


Both  formulas  can  be  written 


~  bp 

a  =  Sa  ,  +  Z  — p  $(t  ,a  .) 
— n  — n-1  —  p!  n'-n-l 


where 


S  =  (I  -  1  Ap)  A  / 


Z  5 


Z  for  corrector  only  methods 


M 

u..£  for  P  (EC)  ‘  methods  , 
M— 


and 


u  -  1+  (1—  £  )  +  ...  +  (1  —  £  ) 
m  P  P 


m-1 
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The  dependence  of  $  on  h  is  not  shown  in  the  notation. 


Definition  of  multivalue  method; 

Let  A  be  a  k+p  by  k+p  matrix  and  Z_  be  a  vector 
of  dimension  k+p.  A  (k+p) -value  method  uses  the  for¬ 
mula 


~  h^ 

a  =  Sa  ,  +  Z  — r  $  (t  ,a  ,  ) 
— n  -n-1  —  pi  n'— n-1 


where  S  is  defined  above  and  $  is  defined  below. 


1. 


Let  0  <  n  <  1.  For  a  corrector  only  method. 


Z  =1,  $(t,a)  =  d>  where 

p  — 


<p  =  F  (t  ,Sa  +  A.  <f>)  ' 


h_P 
P 


and  h  satisfies 


p-q 


X  |r  UilLi 

i=0  p- 


.hP"1  <  i-n  . 


M 

2.  Let  M  be  a  positive  integer.  For  a  P (EC) 

method,  if  £  =  1  then  ^  0  for  some  i  <  p-q,  and 

$(t,a)  =  <j>  (M)  where 


(o)  =  h~PP!  (A^p  ' 


u 


m 


m-n 


^  (m)  _  (1  V*10)\HIWpl  F  (t '  — (j-1) 5 


-  hp 

a  ,  «  =  Sa  +  Z  — r-  4>  ,  x 
-(m)  -  -  p!  (m) 


□ 


- 


M 

P (EC)  methods  are  not  permitted  to  have  £  =  1  and 

P 

£^  =  0  for  i  <  p-q  because  an  explicit  predictor- 
corrector  method  is  really  a  corrector  only  method  if 

£p  =  i. 

Note  that  if  0  <  £  <2  then 

P 


£  =  u  £ 
—  M— 


as  M  -*  °° 


Otherwise  uM  diverges  as  M  +  cof  which  indicates  that  an 
unlimited  number  of  corrector  iterations  can  only  be 
allowed  if  0  <  £^  <  2. 

In  the  case  that  £  =1,  there  is  considerable 

P 

simplif cation  in  the  above  definition: 

<f>  (o)  =  h~PPl  (Aa)  f 

+  <m)  -  F«t'Sa  +  • 

It  has  been  noted  that  using  polynomials  of 
degree  k+p-1  limits  the  accuracy  of  y(t;h).  Thus  it 
is  reasonable  to  construct  y(t;h)  by  means  of  inter¬ 
polation  schemes  which  use  several  of  the  an's  to  define 
y(t;h)  at  any  one  point.  In  order  to  interpolate,  it 

is  necessary  to  assign  a  meaning  to  each  component  of 

c 

a^.  Let  the  correct  value  of  a.n  be  given  by  a  (tn) 
c 

where  a  (t)  is  defined  in  terms  of  the  solution  y(t) 


and  the  stepsize  h.  If  interpolation  schemes  of 


' 
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arbitrarily  high  order  are  used  then  the  accuracy  of 
the  approximate  solution  is  limited  only  by  the  accuracy 
of  the  computed  values  a^ ,  a^,...,^  (and  by  the  number 
of  points  N+l).  Convergence  should  be  redefined  to 
require  that  the  computed  values  converge  in  some 
sense  to  the  correct  values  ac(tn)  so  that  the  defini¬ 
tion  is  independent  of  the  techniques  used  to  construct 
y(t;h).  With  such  a  definition,  there  is  no  a  priori 
limit  on  the  order  of  convergence  of  a  method. 

In  the  case  of  a  k-step  Adams  method,  the  poly¬ 
nomial  p  (t)  is  determined  by  the  components  of 

y  =  col  (y  ,hy' ,hy'  hy '  ,  ,,)  , 

and  thus  there  exists  some  nonsingular  matrix  T  such 
that  a^=  TYn*  Therefore  we  might  define  aC (t)  =  Ty(t) 
where 

y (t)  =  col (y (t)  ,hy  1  (t) ,hy 1  (t-h) , . . . ,hy '  (t- [k-1] h) ) . 


For  k  =  2 ,  we  would  have 


y  (t) 


(1.4)  ac(t) 


hy' (t) 

|  (y' (t)-y' (t-h)) 


A  second  possibility  is  to  define 


aC(t) 


as  the 


vector  of  scaled  derivatives  of  y(t).  For  p  =  1  and 


' 


and  k  =  2  we  would  have 


ac(t)  = 


y  (t) 

hy' (t) 
,  2 


V  y"  <t> 

This  differs  from  (1.4)  by  a  term  of  order  h3. 

For  both  of  the  above  examples  we  can  write 


n 

a  (t)  =  Ay  (t)  where  A  =  col  (A  , A,  , ...  ,A. 


)  is  a 


01'***'  k+p-1 
vector  of  linear  operators.  If  y(t)  is  analytic  then 

each  A^  can  be  expressed  as  a  power  series  in  the 

scaled  differential  operator  hD.  For  example  the 

translation  operator  can  be  written  as  e^D  because 
|*^ 

e  y(t)  =  y (t+h) .  Since  A  indicates  how  to  interpret 
the  ' s ,  we  define  an  interpretation  to  be  a  vector  of 
linear  operators,  each  component  a  power  series  in  hD. 

From  this  point  on,  let  us  adopt  the  Nordsieck 
interpretation 

Ay  (t )  =  col  (y  (t)  ,hy 1  (t)  , .  .  .  ^Y+p-1)  1  y  (k+P~1)  (t) ) 
and  define  a(t)  by 


ai (t)  = 


h^y ^  (t)/i!  for  i  <  r+p-1 


0 


otherwise 


■ 
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Note  that  a(t)  depends  implicitly  on  r,  which  can  be 
adjusted  in  case  some  higher  derivatives  do  not  exist. 

Before  defining  q-convergence ,  let  us  select 
some  norm  | | . | |  for  column  vectors  of  dimension  k+p. 
Corresponding  norms  for  matrices  and  row  vectors  can 
be  defined  in  terms  of  | | . | | .  Define  the  norm 

I  I  •  I  I  H  =  I  I15"1-  I  I  where 


r  i 


H 


0 


hp-q 


o 


h 


p-q 


Definition  of  rth  order  q-convergence  for  multivalue 
methods : 

Let  r  be  a  positive  integer.  A  multivalue  method 
is  q-convergent  of  order  r  if  for  any  Lipschitz  problem 
with  solution  y (t) eCr+^ [0 ,1] , 

| | a  -a (t  ) | |  =0 (hr ) 

i  i  —n  —  n  '  1  jj 

whenever  a^  is  a  function  of  h  satisfying 

E0  =  a.(°)  +  Ofh^P"1) 


□ 


It  is  not  difficult  to  show  that  this  definition  is 
equivalent  to  the  corresponding  definition  for  poly¬ 
nomial  methods  when  r  <  k+1 . 

Section  1.4  Determining  the.  Order  of  a  Method 

In  the  previous  section  the  concept  of  inter¬ 
pretation  was  introduced  in  order  to  show  that  there 
is  nothing  special  about  the  Nordsieck  interpretation. 

It  is  clear  that  the  order  of  convergence  depends  on 
the  choice  of  A,  and  this  suggests  that  for  any  given 
method  we  could  try  to  find  an  interpretation  A*  that 
would  make  the  order  of  convergence  as  high  as  possible. 
The  search  for  such  an  interpretation  uncovers  a  cons¬ 
tructive  definition  for  the  order  of  a  multivalue  method 
Because  y^  ^  appears  in  the  differential  equation,  it 
can  be  shown  that  the  highest  order  attainable  is  k+q+1. 
Therefore  let  us  consider  only  equations  of  the  form 
y  ^  =  f  (t) .  The  following  theorem  asserts  that  for 
any  multivalue  method  satisfying  certain  mild  conditions 
there  exists  a  unique  interpretation  A*  which  makes  the 
method  exact  whenever  f (t)  is  well-behaved  and  h  is 
small  enough. 

The  numerical  solution  of  y^^  =  f  (t)  satisfies 


' 


22 


Given  the  solution  y(t),  the  problem  is  to  find  a* * (t) = 
A*y  (t)  so  that 

a*  (t)  =  Sa*  (t-h)  +  l  y  (p)  (t)  . 

Write  this  as  the  operator  equation 

(1.6)  A*  =  Se-hDA*  +  l  (hD]P 

-  p! 


Theorem  1.1; 

Consider  a  method  for  which  S  has  no  generalized 
eigenvectors  of  rank  >p  corresponding  to  1.  Then  (1.6) 
has  a  unique  solution 


A* 


I  i 

j=0  3 


(hP) 

jl 


/ 


•  oo 

called  the  optimal  interpretation.  For  y(t)eC  [0,1], 
the  power  series  defining  a*  (t)  =  A*y(t)  converges  when 
h  <  ht  where 

ht  =  p(lij^up|y<j)  (t)  |1/j)_1 

and 


p  =  min{  j  log  E,  | :  E,  ^  1  is  an  eigenvalue  of  S}  . 

Proof ;  Substitute  the  complex  variable  z  for  the  opera¬ 
tor  hD  in  (1.6)  to  get 

SA^(z)  =  eZA_(z )  -  py  ez£_ 


(1.7) 


where  A(z)  has  been  substituted  for  A*.  Let 


Qrn(z)  -  qozin  +  q1zII,"1+...+ 


be  a  minimal  polynomial  for  S  and  then  define 

•  *  -i 

Vz)  =  qoz  +  qiz3_  +---+  qj  • 


Repeated  applications  of  (1.7)  yield 


sjx<z>  =  e^XCz)  -|1  i 

p!  4=1 


whence 


Qm(S)  A  (z) 


p  m 

Qffi(ez,X(z)-|r  I  elzQ  (S>I 

3=1  J 


The  left-hand  side  vanishes  leaving 


p  m 

Mz)  = - l  e3Z0.(S  )l  . 

p!Qm(e  )  j=l  3 


Because  S  has  no 
corresponding  to 
of  order  >p,  and 
larities  of  A (z) 


generalized  eigenvectors  of  rank  >p 
z 

1/  Qm(e  )  does  not  have  a  zero  at  z =  0 
so  A_(z)  is  analytic  at  z=  0.  The  singu- 
are  in  the  set 


{log  £  : 1  is  an  eigenvalue  of  S}  . 


Therefore  A_(z)  has  a  MacLaurin  series  with  radius  of 
convergence  p  r  □ 


Consider  for  example  Euler's  method: 


The  optimal  interpretation  is 


f  1  ] 

r  1  i 

r  1  i 

r  1 ) 

J. 

2 

6 

9 

30 

+ 

hD  + 

(hD)  Z 

2 

0 

V  ^ 

1 

0 

0 

A*e 


t 


is  always  easy  to  find,  and  for  this  example 


f  /  -i  — h  v  —  1  ,  > 

(1  -  e  )  h 

h 


In  Section  2.2  the  order  of  a  method  is  defined 

to  be  the  extent  to  which  the  optimal  interpretation 

agrees  with  the  Nordsieck  interpretation.  For  q  >  1 

it  is  necessary  to  impose  a  second  condition  to  ensure 

that  the  difference  between  the  numerical  solutions 

with  starting  values  a(0)  and  a*  (0)  respectively  does 

M 

not  grow  as  n  increases.  Also  for  P(EC)  methods,  it 
is  necessary  to  require  that  the  predictor  matrix  A  be 
sufficiently  close  to  the  Pascal  triangle  matrix  A. 
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) 


Section  1.5  The  Relationship  to  Multistep  Methods 

Let  us  modify  the  proof  of  Theorem  1.1  in  order 
to  show  that  the  numerical  solution  of  a  corrector 
only  method  satisfies  a  difference  equation  of  the 
type  that  is  used  by  multistep  methods.  One  step  of 
a  corrector  only  method  is  defined  by 


a  =  Sa  +  Z  — r  f 
— n  — n-1  —  pi  n 


where  fR  =  f  (t„  ,y„  ,y^ , .  .  .  ,y/P~q) 


n,u  n'“n- 


n 


) .  Thus 


h^ 

Sa  =  a  ,  -  Z  — r  f  .. 
—v  — v+1  —  p!  v+1 


and 


i-i.  hP 


'a  =  a  .  -  Y  S: 

— n-m  — n-m+i  — 


'm+D  ^=i  “  P*  n-m+Z 


Let  Qq (z) ,Q^ (z) , . . . ,Qm (z )  be  the  polynomials  defined 
in  the  proof  of  Theorem  1.1.  Then 


0 


(S)a 


■n-m 


m 


m-1 


Qj  (s>  A 


f„-3  * 


Since  the  pth  row  of  S  is  all  zero,  S  must  be  singular, 
and  so  q  =  0 .  Therefore 


m-1  ,i  ...  m-1  ,p 

(1.8)  T  q.  rr  y  1  •  =  y  (Q.(S)A).  f 

jt0  VD  ~  i  Pi  n-3 


for  i  =  0,1,..., p-q.  What  has  been  derived  in  (1.8)  is 


the  formula  for  an  (m-1) -step  method  where  m  <  h+p. 
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Any  polynomial  Q(z)  ^  0  for  which  Q(S)  =  0  could 
be  used  instead  of  Qm(z),  and  the  numerical  solution 
would  satisfy  the  associated  difference  equation.  Of 
all  difference  equations  derived  in  this  manner,  (1.8) 
has  the  lowest  order. 

Let  T  be  any  (k+p)x(k+p)  nonsingular  matrix  with 
T  T 

6_^T  =  for  i  =  0 , 1 ,  .  .  .  ,p-q  ,p ,  and  consider  the  multi¬ 
value  method  determined  by  A'  =  TAT  ^  and  V0  =  T£_. 

Then  the  associated  multistep  method  given  by  (1.8) 
remains  unchanged  even  though  the  order  of  consistency 
of  the  two  multivalue  methods  may  differ. 

Gear  (1971:151)  states  that  "the  important  charac¬ 
terization  of  a  method  is  the  number  of  values  that  are 
saved  from  step  to  step  and  not  the  number  of  steps." 
Perhaps  this  should  be  amended  to  say  that  the  important 
characterization  of  a  method  is  the  number  of  "independent" 
values  that  are  saved.  We  have  shown  that  a  (k+p) -value 
formula  in  normal  form  is  equivalent  to  a  2 (k+p-1) -value 
multistep  formula.  In  both  cases  there  are  k+p-1 
"independent"  values.  The  "dependent"  values  are 
related  to  these  through  the  differential  equation. 

The  importance  of  the  number  of  "independent"  variables 
is  introduced  by  Gear  (1971:141)  in  connection  with 
stability . 


CHAPTER  2 


NECESSITY  OF  THE  CONDITIONS  FOR  q-CONVERGENCE 

IT  - 1 

A  method  is  said  to  be  q-convergent  like  o (h  ) 
if  it  satisfies  the  definition  of  rth  order  q-conver- 

IT  IT  —  1  • 

gence  with  0 (h  )  replaced  by  o (h  ) ;  that  is , 

I  la  -a(t  )  I  I  =  o (hr  whenever  aA  =  a(0)  +  0 (hr+^  1) . 
In  the  first  three  sections  of  this  chapter  it  is  shown 

r — i 

that  a  method  which  is  q-convergent  like  o (h  )  must 
necessarily  satisfy 

1.  conditions  on  the  eigenvalues  of  the  matrix  S. 

2 .  conditions  on  the  accuracy  of  the  method  when 
f  =  f  (t)  . 

3.  additional  conditions  on  the  accuracy  of  the 
method  when  f  =  f(t,y,...,y^  ^  )  . 

In  the  fourth  section  these  conditions  are  applied  to 
study  multivalue  methods  of  maximum  order. 

In  Chapter  3  these  necessary  conditions  are  shown 
to  be  sufficient  for  q-convergence  of  order  r,  and  in 
Chapter  5  this  result  is  extended  to  some  variable  step- 
size  methods.  For  Chapters  2,  3,  and  5  the  following 
notation  is  very  useful: 
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E  is 


6  .  .  , 
ID 


6  .  is 
-D 


T  . 
is 


the 


(k+p) x  (r+p) 


matrix  defined  by  E. 

ID 


the  jth  column  of  E, 
the  ith  row  of  E, 


P  is  the  (r+p) x (r+p)  Pascal  triangle  matrix 
defined  by 

D 
i 

Notice  that  a(t)  =  EY (t)  where 


Y  (t)  e  col  (y (t) ,hy '  (t) , . . . ,hr+p~1y ^r+p  ^  (t)/  (r+p-1) I ) . 


Also,  if  y (t) eCr+P [0 ,1]  then  by  employing  a  truncated 
Taylor's  series  about  t-h  to  approximate  each  component 
of  Y ( t ) ,  it  follows  that 

Y  (t)  =  PY  (t-h)  +  0(hr+P) 


Section  2.1  The  q-root  Condition 

The  stability  properties  of  a  multivalue  method 
depend  almost  entirely  on  the  matrix  S.  After  all  the 
necessary  conditions  for  q-convergence  are  proved,  it 
is  shown  that 


where  Sq  is  a  (p-q)x(p-q)  upper  triangular  matrix  with 

l's  on  the  diagonal  and  is  a  (k+q)x(k+q)  matrix 

d 

whose  elementary  divisors  have  the  form  (z-£)a  where 
either  |£|  <  1  or  |£|  =1  with  d  £  q.  In  this  section 
we  prove  that  if  the  method  is  q-convergent  then  S 
satisfies  the  q-root  condition  as  given  by  the 

Definition  of  q-root  Condition: 

A  multivalue  method  satisfies  the  q-root  condi¬ 
tion  if  the  elementary  divisors  of  S  have  the  form 
(z-£)^  where  |£j  £  1  and  if  |^|  =1  then  d  <  p.  Further¬ 
more,  if  there  is  an  elementary  divisor  (z-£^)P  where 

d 

=1  then  all  other  elementary  divisors  (z-£) 
with  |^|  =1  satisfy  d  £  q.  □ 

The  q-root  condition,  discovered  by  Gear,  is 
very  similar  to  the  root  condition  of  multistep  theory. 

In  fact,  the  q-root  condition  implies  that  the  associa¬ 
ted  difference  equation  (1.8)  satisfies  the  root 
condition,  which  means  that  the  numerical  solution  of 
a  convergent  corrector  only  method  is  also  the  solution 
of  a  stable  multistep  method. 

Theorem  2,1: 

The  q-root  condition  is  necessary  for  q-conver- 
gence  like  o  (hr  ^) . 


Proof :  Consider  the  most  trivial  Lipschitz  problem 

y(p)  =0  with  yQ  =  Yq  =  ...  =  y^p_1)  =  0.  Let  v  be 
an  arbitrary  vector  and  set  a^  =  hr+p_1v.  The  computed 
solution  =  hr+p  ^Snv.  Hence  q-convergence  like 

3T*™  1 

o  (h  )  implies  that 
(2.1)  || hpSnv|  | H  — ►  0  as  h  — »  0 

The  q-root  condition  is  proved  in  three  parts. 

Part  I :  The  eigenvalues  of  S  are  on  or  inside  the  unit 
circle.  Assume  this  is  false.  Then  there  exists  an 
eigenvector  x  corresponding  to  an  eigenvalue  £  where 
|£|  >1.  Hence  SNx  =  £Nx  and  so  |  |hpSNxj  |  -►  «  as 

h-*  0,  which  contradicts  (2.1). 


Part  II:  An  eigenvalue  of  S  on  the  unit  circle  may  not 

correspond  to  an  elementary  divisor  of  degree  greater 

than  p.  Assume  this  is  false.  Then  there  exists  a 

generalized  eigenvector  x  of  rank  p+1  corresponding  to 

P 

an  eigenvalue  E,  where  \E,\  =1.  Defining  x.=  (S-£I)  p~:Lx 
allows  us  to  write 


-P 


,N 


-P 


fN 

n' 

1 — 1 
l 

Z 

j 

n' 

0 

£  x  + 
-P 

1 

£  X  ,+...+ 

-P-1 

.P. 

rN-p 

s  % 


whence 


' 
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h?sN~  -  £ 


N-p 


^p  p!  ^0  ’  - 


x.  +  O (h) 


Therefore  I  I  h^*SNx  I  I  does  not  tend  to  zero  as  h  -*  0 , 

i  i  — p  H  ' 

which  contradicts  (2.1). 


Part  III;  If  there  is  an  elementary  divisor  (z-£^) 

where  | |  =  1  then  all  other  elementary  divisors 
d 


(z -£)  with  |£|  =1  satisfy  d  <  q.  Assume  this  is 
false.  Let  be  a  generalized  eigenvector  of  rank 

p  corresponding  to  and  let  w^  be  a  generalized 

eigenvector  of  rank  q+1  corresponding  to  an  eigenvalue 


£  on  the  unit  circle  such  that  xA,x, ,...,x  ,  ,  w^ .w, , 

—0—1  — p-1  —0—1 

p-l-i 


.  .  ,w  are  linearly  independent.  Here  x^e(s-^I) 


X  T 
-P-1 


and  wis(S-5Dq-1wq.  For  an  arbitrary  vector  v,  (2.1) 
implies  that 


(2.2)  (SNv)  ±  =  o(Nq) 


if 


i  >  p-q 


Clearly 


SNx  n 

-p-1 


N 

p-1 


:N-p-l 


2o  + 


+ 


N  N-q 

cqj^i  -p-q  '  - 


x  +  o  (Nq) 


If  this  is  substituted  into  (2.2)  with  v  =  x  then 

it  is  not  difficult  to  see  that 


(£0>i  =  <*!>!  =  ••• 


=  (x  )  .  =  0 

-p-q-ii 


for  i  >  p-q 


Similarly 


SNw  = 

-q 


N 

iq 


^N’"qw 0  +  o(Nq) 


and  after  substituting  this  into  (2.2)  with  v  =  w^, 
we  conclude  that 


(Wo’i  =  0 


for 


i  £  p-q 


Therefore  each  of  xA,...,x  , ,wA  is  a  linear  combina- 

—0  — p-q-1  — 0 

tion  of  6  A  ,  6-,  , .  .  .  ,  _  .  .  which  is  impossible  because 

— U  “1.  — p-q-l 


x 


±0' *  *  *  * —p-q 


f-->-  -l'^O  are  linearlY  independent 


□ 


Section  2.2  The  Accuracy  of  the  Corrector 

In  this  section  two  conditions  are  derived  which 
are  necessary  for  convergence  when  f  =  f  (t) :  the  cor¬ 
rector  condition  and  the  growth  condition.  These  two 
conditions  determine  the  order  of  the  linear  part  of 
the  formula,  Sa  ,  . 

In  the  next  section  the  predictor  condition  is 
derived,  which  is  necessary  for  convergence  when 
f  =  f  (t,y , . . . ,y ) .  The  predictor  condition  deter¬ 


mines  the  order  of  the  nonlinear  part  of  the  formula, 
hi 
P 


~  V-i  P 

i  pT  $(tn'^n-lK 


These  conditions  are  combined  in  the 


■ 


Definition  of  rth  q-order  condition: 


Let  r  be  a  positive  integer.  A  multivalue  method 
satisfies  the  rth  q-order  condition  if  S  has  no  genera¬ 
lized  eigenvector  of  rank  >p  corresponding  to  1  and  if 
the  method  satisfies 

1.  the  corrector  condition: 


Dij  =  5ij  for  j  <  min (i,p-q)+r 

where  D  is  the  (k+p)x(r+p)  matrix  defined  by 

D^j=  (dj )  ^  or  equivalently  D  =  row  (d^d^,..., 

d  ,  ,  )  and  the  d.'s  are  defined  by  Theorem  1.1. 

— r+p-1  — j  J 

2.  the  growth  condition:  for  r+p-q  £  j  <  r+p-1  and 

i  >  p-q* 


3. 


(Sn  (d  .  -  6.))  .  =  0(nj-r  P+q)  . 

— j  — j  i 

M 

the  predictor  condition:  for  a  P  (EC)  method  if 

£  =1  then 

P 


for  j  <  r+p-Mc 


where  c  is  the  smallest  integer  >q  such  that 

Z  f  0  . 

p-c 


□ 


It  is  easy  to  see  that  the  (r+l)th  condition  implies  the 
rth  q-order  condition.  The  corrector  condition  states 


that  D  has  the  form 


t 

1 


1 


col 

col 

col 

r 

r+p-q 

r+p-1 

x 

.  .  .  X 

.  .  X 

• 

• 

• 

• 

• 

•  • 

• 

• 

0 

X  X 

.  .  x 

c 

X 

e  .  X 

row 


p-q 


0 


1 


x 


x 

/ 


This  implies  that  |  | a*  (t) -a  (t)  |  | R  =  0  (hr)  for  small  h. 

For  q  =  1  the  growth  condition  does  not  apply;  other¬ 
wise  it  ensures  that  |  |  Sn  (a  (0)  -a*  (0))|  |  =  0(hr). 

The  necessity  of  parts  1  and  2  of  the  rth  q-order 
condition  is  proved  by  considering  the  problem  y  ^  =  e ^ 
with  y^  =  X1  where  |x|  <  1.  Then  a*  (t)  exists  for  h  <  p. 
Since  a*  (tg )  ,a*  (t-^)  , .  .  .  ,a*  (tN)  is  a  numerical  solution. 


' 

it  follows  that  for  any  other  numerical  solution 


(2.3)  an_2*(tn)  =  Sn(aQ-a*(0))  . 


Theorem  2.2: 

Parts  1  and  2  of  the  rth  q-order  condition  are 
necessary  for  q-convergence  like  o  (h  ). 


Proof:  By  definition. 


CO 


a*  (t)  =  l  d.  hV]’  (t) 

j=o  ~3  31 


(2.4) 


=  DY (t)  +  o (hr+P  1) 


Let  the  starting  value  be  correct:  an  =  a(0)  =  EY(0). 

— u  —  — 

r — i 

q-convergence  like  o (h  )  means  that 


a^  =  a (1)  +  Ho (hr-1) 

(2.5)  =  EY  (1 )  +  Ho  (hr-1)  . 

Substituting  (2.4)  and  (2.5)  into  (2.3)  with  n  =  N  gives 
EY (1 )  +  Ho (hr_1)  -  DY (1 )  -  o  (hr+P_1) 

=  SN(EY(0)  -  DY  (0 )  -  o(hr+P~1))  . 


36 


It  follows  from  (2.1)  in  the  proof  of  Theorem  2.1  that 
hpS*  *v  -  Ho(l)  for  any  arbitrary  vector  v,  and  hence 
SNo (hr+p  "*")  =  Ho  (hr  1)  .  Therefore 

SN(D-E)Y(0)-  (D-E)Y(l)  =  Ho  (hr_1)  , 

or  equivalently 


r+jD-i 

j=0 


XjSN  (d.-6  .)  -  ?  AjeX  hD  r  - - /’-r“1 


-3  ~3 


j=0 


,  .  (d.-6.)  =  Ho  (h  -1-)  . 

J  *  3  3 


The  functions  l,A,...,Ar+p  1,  eX,  AeX,...,Ar+p  1eX  are 
linearly  independent  on  [-1,1],  and  so  there  exists  a 
set  J  of  2 (r+p)  points  in  [-1,1]  such  that  these  func¬ 
tions  are  linearly  independent  on  J.  Therefore  each 
of  the  2  (r+p)  coefficients  of  these  functions  must  equal 
Ho  (hr  ■*■)  .  Thus 


D4J  =  6,,  +  o(hr-l-j  +min(i,p-q)) 


ID 


ID 


from  which  the  corrector  condition  follows .  Also  for 

i  >  p-q 

<SN(d.-S.)).  =  o(N^r+1-P+<3) 

— D  — D  i 

where  N~X  has  been  substituted  for  h.  Since  the  eigen¬ 
values  of  S  are  on  or  inside  the  unit  circle,  either  the 
left-hand  side  converges  to  zero  as  N  00  or  there  is  some 


I  .  I 
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nonnegative  integer  m  such  that  N  m  (SN  (d_. -6__. )  )  ^  is  a 
bounded  nonconvergent  sequence.  From  this  it  follows  that 
the  right-hand  side  can  be  replaced  by  0  (N-1  r  P+c*)  if 
j  £  r+p-q  thus  proving  the  growth  condition.  □ 


Section  2.3  The  Accuracy  of  the  Predictor 

M 

This  section  applies  only  to  P (EC)  methods. 

M 

Recall  that  the  definition  of  a  P  (EC)  method  with  £  =  1 

ir 

required  that  ^  0  for  some  i  1  p-q.  This  restriction 

was  made  to  ensure  the  existence  of  c,  which  is  the 

smallest  integer  >  q  such  that  £  ?  0.  If  c>q  then  the 

p  c 

method  is  partially  explicit  because  the  predicted  value 

of  y  'P"^)  j_s  left  unchanged  by  subsequent  corrections. 

The  predictor  condition  requires  that  y^  (o) ' 

y'  / . . . /Y  have  order  of  accuracy  r  -  Me ,  and  it 

Jn,  (0)  (0) 

suggests  that  each  correction  increases  the  order  by  c. 

If  £p  ¥  If  the  predicted  values  have  order  r-q  and  this 
is  not  increased  until  the  last  correction,  when  it  is 
increased  by  q.  Gear  (1971:191)  states  in  one  of  his 
proofs  that  if  the  corrector  has  order  r  for  1-differen¬ 
tial  equations  then  the  predictor  must  have  order  r-M. 
This  statement  is  not  true  if  £p  =  1  and  £p-1  =  °- 

To  prove  the  necessity  of  the  predictor  condition, 
we  consider  the  problem  y^^=  y^  with  Yq  =  •••  = 


,  (P-D  = 
0 


1. 


' 


Theorem  2.3: 


Part  3  of  the  rth  q-order  condition  is  necessary 

r — i 

for  q-convergence  like  o (h  ) . 

Proof:  Assume  £  =  1  and  r+p-Mc >  0;  otherwise  there  is 

-  P 

nothing  to  prove.  The  major  thrust  of  the  proof  con¬ 
sists  of  arriving  at  a  contradiction  by  assuming  that 
the  predictor  condition  is  violated.  That  is,  a  con¬ 
tradiction  is  obtained  by  assuming  that  an  integer  s  <  r 
exists  such  that 

(2.6)  A  .  =  P  .  for  j  <  s+p-Mc 

PD  P3  J 

(2.7)  Ap, s+p-Mc  =  Pp, s+p-Mc  +0t 

where  a  ^  0.  Let  us  consider  the  following  two  initial 
value  problems : 

1.  y^  =  et  with  yQ =  . . . =  y^p  =  1 

2.  y^  =  y  (p-c)  with  yQ=  .  . .  =  y^p  =  1  , 

which  have  the  same  solution  y(t)  =  efc.  By  the  hypo¬ 
thesis  the  numerical  solution  of  problem  1  satisfies 


and  the  numerical  solution  of  problem  2  satisfies 


(2 . 8 )  an  =  EY(t  )  +  Ho  (hr-1 


Subtracting  these  equations  shows  that  the  difference 
c  E  a  -a'  satisfies 


— n  — n 


r-1 


c  =  Ho  (h 


r 


and  in  particular 


(Sn>  o  =  °<hr  1>  • 


The  remaining  part  of  the  proof  consists  of  finding  re¬ 
currence  relations  that  are  satisfied  by  the  numerical 
solutions  an  and  a/.  Then  by  using  assumptions  (2.6) 
and  (2.7),  it  is  shown  that  the  difference  ( ) q  satis¬ 
fies 


with  a  ^  0  and  &p_c  ^  0-  The  two  equations  for  (c^q 
are  contradictory  because 


lim 

h+0 


i=0 


' 


. 


■  • 


Thus  assumptions  (2.6)  and  (2.7)  are  false  and  the 
method  satisfies  the  predictor  condition.  To  complete 
the  details  of  the  proof  it  must  be  shown  that 

(2.10)  a'  =  SaT  +  i  ^  enh 

— n  — n  —  pi 

(2.11)  a  =  Sa +  £  ~  (enh  +  hS(3enh  +  o(hs)) 

1 1  11  M  • 

and  thus 

£n  =  s£n_i  +  i  |t  (hSBenh  +  o(hs))  . 

This  recurrence  for  £n  yields 

o  =  hs  ?  Sn-i  J  ^  pe5h  +  y  o(hs+P)Sn-3jl  . 
j=l  ~  P!  jil 

Once  equations  (2.10)  and  (2.11)  have  been  verified,  it 
is  easily  checked  that  the  first  summation  is  the  nth 
term  of  the  numerical  solution  of  the  problem 

6(p)(t)  =  ^  with  6Q  =  ...  =  6qP_1)=  0  . 
Since  the  method  is  q-convergent 


‘^>0 


=  h‘ 


<6(tn) 


+  o(l))  + 


n 

l 

j=i 


(hs+P) 


(Sn  JZ) 


Also,  by  Theorem  2.1  the  method  satisfies  the  q-root 
condition  and  thus  (Sn  ^) q  =  0 (h^  p) .  Therefore 
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(£n)0  =  h  6(tn)  +  o(h  }  ' 


which  yields  (2.9).  Equation  (2,10)  may  be  verified 
by  using  (1.5).  This  leaves  us  to  do  the  Verification 
of  (2.11):  By  definition 


hP 

(2.12)  a  =  Sa  ,  +  £  — r  cj> 

— n  — n-1  —  p!  Yn,  (M) 


where 


(2.13)  <p 


n,  (0) 


=  h  P  p!  6TAa  ,  , 

c  — p  — n-1  ' 


(2.14)  <j>  ,  ,  =  h  P+C(p- 

yn, (m)  ^  -p 


T  hp 

(p-c)  !  6  (Sa  ,  +  £  ~  cf)  ,  ,  .) 

c  -p-c  -n-1  —  p!  Yn, (m-1) 


Using  £^=  1,  it  follows  that  (I-  £_  _6p)£_  ”  0'  and  thus 
T  M  T 

(I  -  £  6  )  =  I  -  £  6  .  Therefore  multiplying  (2.12)  by 

—  — p  —  — p 

I  -  £_  6_p  gives 


Sa  ,  =  (I  -  £  6T)a 
— n-1  - p  — n 


hP  nh 


r-1 


=  EY(tn)  -  £  e““  +  Ho  (h  by  (2.8) 


Substitute  this  into  (2.14)  to  get 


/o  ncr\  ,  nh  ,  c  (p-c)!  .  , ,  nh.  ,  r-1. 

(2.15)  cb  .  v-  e  =  h  — — ; £  (6  ,  1  .  -e  ) +o  (h  ) 

Yn, (m)  p!  p-c  Yn, (m-1) 


From  the  initial  assumptions  (2.6)  and  (2.7) 


(6TAE) .  =  ( (eT  +  aeT  )P). 

— p  j  — p  —s+p-Mc  3 


for  j  <  s+p-Mc 


whence 


6^AEY(t  )  =  ^  enh  +  oih(  -J  enh  +  0  (hS+P“Mc+1 ) 

— p  —  n-1  p!  (s+p-Mc) !  ' 


^s+p-Mc 


Therefore  from  (2.8)  and  (2.13) 


+n,{0)  =  h-Pp!^A(EY(tn_1)  +  Ho  (h^1)  ) 


,  ,  S-Mc  .  , 

nh  .  _  h  p !  nh 


=  e  +  a 


(s+p-Mc) i 


****  ,  ~  ,,  s-Mc+1. 
e  +  0  (h  ) 


+  o  (h“p+rnin (p-q/ s+p-Mc )+r-l} 


and  so 


s— Me 

,  nh  h  ap!  nh  ,  ,,  s-Mcx 

(2.16)  4>  rn\  ~  e  =  i — : - 77  \  ,  e  +  o  (h  ) 

rn,  (0)  (s+p-Mc) ! 


n  H 

Solve  for  <J>n  ^  -  e  using  (2.15)  and  (2.16) 


,  nh  ,  sfl  nh  .  sx 

-  e  =  h  6e  +  °(h  )  * 


The  proof  of  (2.11)  follows  from  this  and  (2.12). 


□ 


Section  2.4  Methods  of  Maximum  Order 

The  following  lemma  partially  characterizes  the 
S  matrix  when  the  rth  q-order  condition  is  satisfied. 
Furthermore  it  is  used  to  gain  some  insight  into  the 
vectors  dn ,d.  , . . . ,d  ,  . 

—0—1  — p-1 


Lemma  2.1: 


If  S  has  no  generalized  eigenvector  of  rank  >p 
corresponding  to  1  then 


(2.17)  SD  =  DP  -  £  e  P  . 

- p 


Proof :  From  (1.7)  in  the  proof  of  Theorem  1.1, 


SA_(z)  =  eZX(z)  -  eZ£ 

ir  • 


with 


00  D 

z  J 


i<z>  =  l.  ij  j. 


j  =  0 


Take  the  jth  derivative  where  0  1  j  £  r+p-1: 


SX  (j) (z)  =  ez  l 


c  •  ^ 
3 

X  (l)  (z)  -  £  ez  l 

f  . 

3 

t- 

o 

II 

•H 

1 

1 

-L 

d  z 


dz 


i  p! 


Set  z  =  0  to  get 


Sd.  = 

~3 


r  ,  > 

3 

d.  - 

r  ,  *v 

3 

=  0 

i 

-3 

IPJ 

£ 


thus  verifying  the  jth  column  of  (2.17) 


Corollary: 

If  a  method  satisfies  the  q-root  condition  and 
the  rth  q-order  condition  then 


S.  .  =  (  ~  .^)A)  .  .  for  j  <  min (i+1 ,p-q) +r  . 

-L  J  lr  J-  J 


Proof :  Rearrange  Lemma  2.1  to  get 


SE  =  (I  -  l  6T)EP  +  (D-E)P  -  S(D-E)  . 

—  — p 

Column  by  column  it  can  be  shown  that  S  has  the  required 
form.  □ 


If  the  hypotheses  of  the  corollary  are  satisfied  with 
r  =  1,  then  S  has  the  form 


Consider  a  method  which  is  q-convergent .  By 
Theorem  2.1  the  q-root  condition  must  hold,  and  so  S 
has  no  generalized  eigenvectors  of  rank  >p  correspond¬ 
ing  to  1.  Hence  by  the  last  equation  in  the  proof  of 
Lemma  2.1, 
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(S-  I)dj 


for  j  <  p-1  , 


or  in  other  words 


(S  e  <dQ ,<3^, . . . 

where  angle  brackets  denote  the  linear  span  of  the 
enclosed  vectors.  From  this  it  follows  that  for  j  = 
0,1,..., p-1,  dj  is  a  generalized  eigenvector  of  rank 
<j+l  corresponding  to  1. 

The  next  lemma  shows  that  it  is  not  necessary 
to  determine  D  in  order  to  verify  the  corrector  condi- 
tion  and  the  growth  condition.  Any  matrix  D  which 
satisfies  (2.19)  can  be  used  in  place  of  D  in  the  rth 
q-order  condition. 

Lemma  2.2: 

Assume  a  method  satisfies  the  q-root  condition, 
and  suppose  that  there  exists  a  (k+p)x(r+p)  matrix  D 
such  that 

(2.18)  6^^  =  for  j  <  min(i,p-q)+r  , 

(2.19)  S8  =  DP  -  l  eT  P  , 

- p 

(Sn(D~E)).j  =  0(n^r-P+C2) 


and 


. 
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for  r+p-q  <  j  <  r+p-1  and  i  >  p-q.  Then  the  method 
satisfies  parts  1  and  2  of  the  rth  q-order  condition. 

Proof :  The  proof  is  divided  into  two  parts.  The  first 

part  establishes  two  important  facts  (equations  (2.20) 
and  (2.21))  about  certain  subspaces  invariant  under  S. 
The  second  part  actually  proves  the  statement  of  the 
lemma. 

Part  I:  Let  X^_^  t>e  the  subspace  of  generalized  eigen¬ 
vectors  of  S  corresponding  to  1  and  define  the  subspaces 


X.  =  (S-I)p  1  1  X  -j 

l  p— 1 


Let  d.  =  De  .  .  We  assert  that 
—J  -D 


(2.20)  cl.  e  X. 

-3  3 


for  j  <  p-1 


This  assertion  is  verified  by  successively  proving  the 
following  statements: 


p-1  — 1  p-i  —2  p-1 ' 


dn  eX 


d»eX 


*  •  •  t 


As 


W  ^i£Xi 


bi 
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In  order  to  prove  that  d . eX  _  ,  use  the  following 

1  p  j. 

information : 

1.  the  ith  column  of  (2.19): 


_ 

f  •  > 

1 

A 

d^  ^  +  ...  + 

r  ,  \ 

l 

i-1 

0 

A 

^0 


2*  -O' *  *  * '-i-leXp-l  ‘ 


In  order  to  prove  that  d . eX .  where  j  <  p-1 ,  use  the  fol- 
lowing  information: 

1.  the  (i+l)th  column  of  (2.19): 


,  _  % 

l+l 

A 

d.  + 

i+1 

d .  +  ...  + 

i+1 

i 

—i 

i-1 

-l-l 

0 

— u 


2. 


d  •  i  eX  .  -| 

—l+l  j+1 


3. 


—0 '  *  *  * i— l^Xj  * 

It  is  clear  from  (2.18)  that  for  j  £  p-q 


d.  =  <5.  +  some  vector  in  <5~,. > 
-3  -J  — 0  —j -r 


By  induction  on  j  it  follows  that  for  j  <  p-q 


<Aq  '£i  '  •  *  *  'A-i  >  -  <^0  '  —  1 '  *  *  *  ■ i> 


c  XAuX, u . . . uX . 

-01  j 


by  (2.20) 


c  X. 

“  3 


<  U 
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Since  ^  {()},  there  is  a  generalized  eigenvector  Xp_i 

of  rank  p  corresponding  to  1.  Let  Xj_  =  (S-I)p  — p-1* 

By  the  q-root  condition,  X  _ ,  is  the  direct  sum  of 

<x^ , •  *  *  *  ' 2£p_]_>  aRd  the  subspace  of  generalized 

eigenvectors  of  rank  <  q.  Thus  dim  X  _  =  dim 

p  q  J. 


(S-I)M  Xp-i  -  P”q*  and  since  dim  X^_^  <  dim  X^  ,  we  have 
dim  Xj  <  j+1  for  j  <  p-q.  Therefore 


(2.21)  X.  =  «5n  6n  , .  . .  ,6  .  > 

3  -0,-1'  -j 


if  j  <  p-q  • 


Part  II:  Subtract  (2.19)  from  (2.17) 


S (D-D)  =  (D-D) P  . 

A 

Let  v.  =  (D-D)  e  .  ,  for  u  =  —  r,— r4"l,...,0,...,p— 1.  Then 
—J  ~J+r 

by  the  same  argument  as  in  (2.20)  it  is  possible  to  show 
that 


v .  eX . 
—J  3 


By  (2.21),  v .  e  <  i0'-l'* 


i 


0 


. ,6^  >  for  j  <  p-q  so  that 
if  i  >  j  and  j  <  p-q  . 


Equivalently 


0  for  j  <  min(i,p-q)+r 


The  corrector  condition  follows  from  this  and  (2.18). 

Let  p-q  <  j  <  p-1.  Since  v_.  eX ^  ,  there  exist  vectors 

x.eX.  such  that 
—1  1 
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Snv .  =  v .  +  n 
-D  -D 


( 


+  some  vector  in  < 


Therefore  when  i  >  p-q 


(Sn(D-E))ij  =  (Sn(D-D))_  +  (SR  (D-E )  )  ± j 


which  proves  the  growth  condition. 


□ 


As  an  example  it  can  be  shown  that  parts  1  and  2 
of  the  rth  q-order  condition  are  satisfied  when 


(  • 


for  j  <  r+p-1 


Simply  apply  Lemma  2.2  with  D  =  E.  Recall  that  the 
derivation  of  the  multivalue  method  starts  with  the 


Pascal  triangle  matrix  A. 

Gear  (1971:203)  shows  that  there  exist  strongly 
stable  (k+p) -value  methods  which  are  q-convergent  of 
order  k+q.  A  strongly  stable  method  has  no  extraneous 
roots  on  the  unit  circle.  In  fact  if  q  =  1,  £  can  be 
chosen  so  that  the  extraneous  roots  are  all  zero.  Gear 
(1971:154)  gives  coefficients  for  some  of  these  higher 
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order  analogues  of  Adams'  method.  The  case  p  =  1  and 

k  =  5  is  Nordsieck's  method. 

There  are  no  (k+p) -value  methods  which  are  q- 

convergent  of  order  k+q+1  when  k+q+1  is  odd  for  the 

cases  p  =  q  =  1  and  p  =  q  =  2.  Let  us  prove  this 

statement  by  assuming  that  such  a  method  exists.  Set 

r  =  k+q+1  in  Lemma  2.1.  If  the  corrector  condition  is 

satisfied,  the  (p,k+p)th  component  of  (2.17)  forces 

M 

&  =  1.  Clearly  the  order  of  a  P  (EC)  method  with  £  =  1 

p  p 

cannot  exceed  the  order  of  the  corresponding  corrector 
only  method,  and  to  every  corrector  only  method  there 
corresponds  a  stable  (k+q-l)-step  method.  However,  for 
the  two  cases  under  consideration  it  is  known  (Henrici 
(1962))  that  the  order  of  a  stable  (k+q-l)-step  method 
cannot  exceed  k+q  if  k+q-1  is  odd. 

There  exist  (k+p) -value  methods  which  are  q- 
convergent  of  order  k+q+1  when  k+q+1  =  4  and  q  =  1 ,  or 
equivalently,  there  exist  (p+2) -value  methods  which  are 
1-convergent  of  order  4.  To  show  this,  we  use  Lemma 
2 . 2  with 


p+2 

\-l 

p+2' 

.  P  . 

.  i  . 

A 


A 


i. 

l 


. 
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f 


\ 


1 


0 


1 


0 


1 


0 


d 


Then  all  but  the  last  column  of  (2.19)  are  automatically 
satisfied.  The  last  column  of  this  equation  can  be 
written  in  the  form 


f 


0  x 


0 


x 


X 

-1 

X 


X 


X 

0 

X 


^ - 

o 

C*T 3  •  « 

\ _ 

f  \ 

X 

• 

1 — I 

1 

a> 

•  •  <r0 

= 

• 

X 

A 

d 

P 

0 

d  ^ 

l  P+1  J 

X 

»  4 

Clearly  d^  =  0 .  It  so  happens  that  the  same  value  of 
dp+^  satisfies  both  the  (p-l)th  and  the  (p+l)th  equation. 

A  A  A  _ 

The  other  components  dp_^ ,d^_2 ' • *  * can  determined 
by  back  substitution.  The  value  of  dQ  is  immaterial. 

This  method  is  weakly  stable  because  the  extraneous 
eigenvalues  are  0  and  -1. 

Assume  some  method  satisfies  the  (k+q+2)th  q- 
order  condition.  With  r  =  k+q+2,  the  (p-l,k+p)th  com¬ 
ponent  of  (2.17)  is 


. 

[k+PU  n 

l  P  J  P-1 


k+p 

J>“1 


and  the  (p-1 ,k+p-l) th  component  is 


'k+p+lj  j 

.  P  J  P-1 


k+p+1 

l  P-1 


Because  these  two  equations  are  contradictory,  no  method 
satisfies  the  (k+q+2)th  q-order  condition. 


CHAPTER  3 


SUFFICIENCY  OF  THE  CONDITIONS  FOR  q-CONVERGENCE 

In  the  first  section  of  this  chapter  q-stability 
and  rth  order  q-consistency  are  defined,  and  it  is  shown 
that  together  these  properties  are  sufficient  to  prove 
that  a  method  is  rth  order  q-convergent .  In  the  second 
section  q-stability  is  proved  from  the  q-root  condition 
and  the  1st  q-order  condition,  and  in  the  third  section 
rth  order  q-consistency  is  proved  from  the  q-root  condi¬ 
tion  and  the  rth  q-order  condition. 

For  Chapters  3  and  4  the  computed  solution  a^ , 
a is  defined  by 


Sa  , +  i  — r  $ (t 
— n-1  —  pi 


N 


r  *  •  •  r 


where  r^  is  the  starting  error  and  r ^ , £2,***'£N  are 
errors  arising  from  finite  precision  arithmetic. 
These  errors  are  lumped  together  as  the  round-off 
error 


/  *  •  •  9 


We  measure  R  with  the  norm 
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1C '3 
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maX  I  I  V  qn“jv.  I  I 

0<n<N 1 1 j£0  &  — j 1 1 H  ' 


which  resembles  the  norm  introduced  by  Spijker  (1971) . 


Section  3.1  q-stability  and  q-consistency 

Roughly  speaking , stability  is  convergence  of  order 
zero;  that  is,  a  method  is  stable  if  I  la  |  |  TT  is 

uniformly  bounded  for  all  h  whenever  h-*-“P  |  |  |  |  is  uni¬ 

formly  bounded  and  r^= .  .  .=r^=0_.  Such  a  definition  is  not 
very  useful  for  proving  theorems  or  analyzing  the  effects 
of  round-off  error.  A  much  more  useful  definition  is 
the 


Definition  of  q-stability: 

A  multivalue  method  is  q-stable  if 


1.  for  any  Lipschitz  problem,  there  exists  a  constant 

k^  such  that 


x 

a  -a 
— n  -n 


< 

H  “ 


XX  X 

where  ,a  . . . ,aN 
no  round-off  error 


is  the  numerical  solution  with 
("x"  is  for  "exact"). 


2. 


there  exists  a  constant  kg  such  that 


H_1  Sn 


<  ks  h1"? 


□ 


* 
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For  a  q-stable  method  it  follows  that 


x 

a  -a 
— n  -n 


H 


-  koksh 


1-p 


N 

I 

n=0 


r 

— n 


which  is  a  more  conventional  definition  for  stability 
(see  Chartres  and  Stepleman  (1972) ) . 

The  usual  definition  of  consistency  (Gear  (1971: 
181))  states  that  a(t)  =  a(t)+0(hr+p)  where 

-  hp 

a(t)  =  Sa(t-h)  +  l_  — -  $ (t ,a (t-h) )  . 

P  • 


In  other  words,  the  local  truncation  error  arising  from 


the  correct  solution  a(t)  must  be  of  order  h 


r+p 


Al¬ 


though,  this  definition  is  satisfactory  for  multistep 

methods,  it  is  too  strong  for  multivalue  methods.  A 

method  which  is  rth  order  q-consistent  in  the  usual  sense 

would  have  to  satisfy  D.  .  =  6.  .  for  j  <min(i,p)+r,  which 

ID  ^D 

is  stronger  than  the  corrector  condition. 

At  the  beginning  of  Chapter  1,  the  convergence 

of  a  method  is  demonstrated  by  finding  an  interpretation 

of  the  computed  values  for  which  the  local  truncation 
2 

error  is  0 (h  )  and  which  differs  from  the  Nordsieck 

interpretation  by  0  (h) .  Therefore  a  method  is  defined 

to  be  q-consistent  of  order  r  if  there  exists  an  inter- 

r+ 1 

pretation  for  which  the  local  truncation  error  is  0  (h  ) 


■ 
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and  which  differs  from  the  Nordsieck  interpretation  by 
0  (h  ) .  This  special  interpretation  is  represented  in 
the  following  definition  by  the  matrix  D  in  the  same 
way  that  E  represents  the  Nordsieck  interpretation  and 
D  represents  the  first  r+p  terms  of  the  optimal  inter¬ 
pretation  . 

Definition  of  rth  order  q-consistency : 

Let  r  be  a  positive  interger.  A  multivalue  method 

✓s 

is  q-consistent  of  order  r  if  there  exists  a  matrix  D 
such  that  for  any  Lipschitz  problem  with  solution 
y  (t) £Cr+^ [ 0 , 1] ,  the  function  aS  (t)  =  DY(t)  satisfies 

1.  |  |aS  (t)-a(t)  |  |H  =  0 (hr ) 

2.  |  |Sn(as  (0)-a(0)  )  |  |R  =  0(hr) 

3.  aS  (t)  =  as  (t)  +  0(hr+^) 

,  p 

where  aS  (t)  =  Sas  (t-h)  +  £_  — y  0  (t,aS  (t-h)  )  .  □ 

p  i 

Part  2  of  the  definition  has  been  introduced  for 

the  case  q>l;  it  is  not  necessary  if  q  =  1 .  Part  2 

states  that  the  special  starting  value  differs  from  the 

correct  starting  value  by  0(hr)  when  measured  by  the 

norm  nmaX^ I  I Sn . I  Part  1  can  be  deduced  from  part  2 

u<n<N 11  1 1 n 


. 


- 


and  is  therefore  redundant,  but  it  is  included  to  make 
the  definition  more  meaningful. 

The  symbol  D  is  used  in  the  definition  in  order 
to  suggest  the  correspondence  between  the  three  parts 
of  q-consistency  and  the  three  equations  of  Lemma  2.2. 

In  Section  3.3  it  is  shown  that  part  3  of  rth 
order  q-consistency  is  equivalent  to 

(3.1)  $  (t , aS  (t-h) )  =  y(p)(t)  +  0(hr) 

when  D  =  D,  assuming  D  exists.  If  $  satisfies  a  Lips- 
chitz  condition  as  in  Lemma  3.2  and  if  part  1  of  q- 
consistency  holds,  then  (3.1)  is  equivalent  to 

$  (t  ,a  (t-h) )  =  y(p)(t)  +  0(hr)  . 

The  following  example  shows  that  there  are  very 
useful  methods  for  which  the  usual  definition  of  consis¬ 
tency  is  inadequate.  The  2-step  Adams-Moulton  method  in 
normal  form  has  the  parameters 


'  1 

1 

% 

1 

5/12 

A  = 

0 

1 

2 

and  Z  = 

1 

0 

1 

> 

.  1/2  . 

The  truncated  optimal  interpretation  gives  us 


* 

. 

' 


• 

aS(t) 


1  0  01/4' 


0  10  0 


Y(t)  . 


001  -3/2 


The  3rd  1-order  condition  is  satisfied  since  the  first 
3  columns  of  D  are  identical  to  those  of  E.  However, 
the  usual  definition  of  3rd  order  consistency  cannot 


be  satisfied  unless  D.  .  =  6.  .  for  j  <  min(i,l)+3,  which 

13  J 


is  not  the  case. 

The  next  example  of  a  PtEC)"*"  method  for  second 
order  equations  demonstrates  that  2-consistency  does 
not  imply  1-consistency.  For  the  parameters 


'111' 

r  0  ' 

A  = 

0  0  1 

and  £  = 

0 

t 

0 

1 — 1 

1 

0 

_ f 

2 

v.  > 

be  shown  that 

'ill' 

'  1 

0  -1/6' 

S  = 

0  0  1 

and  D  = 

0 

1/2  0 

0  10 

> 

0 

1/2  1  , 

The  1-root  condition  and  the  1st  2-order  condition  are 
satisfied  but  the  1st  1-order  condition  is  not.  Note 
that  the  definition  of  rth  order  q-consistency  depends 
on  q  through  the  use  of  the  norm  | | . | L,. 


. 


. 


' 


. 
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The  following  theorem,  from  Chartres  and  Stepleman 
(1972),  enables  us  to  establish  rth  order  q-convergence 
by  proving  q-stability  and  rth  order  q-consistency . 


Theorem  3,1: 

Consider  a  method  which  is  q-stable  and  q-consis- 
tent  of  order  r.  Then  the  computed  solution  for  a 
Lipschitz  problem  with  solution  y (t) eCr+^ [0 , 1]  satisfies 

I  I  lH  *  k0l?.ls  +  0(hr)  ' 

and  therefore  the  method  is  q-convergent  of  order  r. 


Proof:  By  the  triangle  inequality 


a  -a  (t  )  I  I „  < 
— n  —  n  1  1  H 


x 

a  -a 
— n  — n 


H 


s  , ,  x  X  |  | 
a  (t  ) -a  L, 
—  n  — n  1  1  H 


a  (t  ) -a  (t  ) 
—  n  — n  n 


H 


where  ar  (t)  satisfies  all  three  parts  of  rth  order  q- 
consistency.  By  q-stability 


a  -a* I  I u  ^  | R | 

— n  — n  1  1  H  0  1  — 1  S 


Notice  that  aS  (tQ ) , as  (t1) , . . . ,as  (tN)  is  the  computed 


solution  with 


. 


and  thus  by  q-stability 


| | a  <tn) “Sn I Ih  S  k0 Is 

£  k0  0<n<NllSn<^S(t0)-^(t0))l 
+  kokshl_P  .^1  l£S<tj)-iS(tj> 

and  by  parts  2  and  3  of  rth  order  q-consistency 

lla^V-a^H  “  ko°<hr)  +  kokSNO(hr+1) 

=  0  (hr )  . 


By  part  1  of  q-consistency 

il?S<V-a(tn)||H  =  0(hr)  . 


□ 


Section  3.2  Proving  q-stability 

The  proof  of  q-stability  requires  that  a  method 
possess  two  important  properties: 

||H""1Sn||  <  ksh1-”p  for  some  constant  kg ,  which 
is  proved  from  the  q-root  condition  and  the  1st 
q-order  condition. 


1. 


■ 


61 


2.  for  any  Lipschitz  problem,  there  exists  a  cons¬ 

tant  L  such  that 

|  $  (t  ,a2) -$  (t,a^)  |  £  L|  l^-a1]  |H  , 
which  is  proved  from  the  1st  q-order  condition. 

The  first  property  is  a  consequence  of  the  following 
lemma,  which  examines  the  size  of  Sn. 

Lemma  3.1: 

If  a  method  satisfies  the  q-root  condition  and  the 
1st  q-order  condition,  then  the  components  of  Sn  are  of 
order 


n 

•  • 

• 

i 

iL 

1 

H 

nP-'S.  . 

p-i 
.  n^ 

• 

p-1 

.  n^ 

• 

•  \ 

« 

• 

• 

1 

• 

C  * 

• 

• 

• 

• 

• 

1 

•  • 

n  ' 

*2 

n 

.  n-3+1. 

• 

q+1 
.  n^ 

0 

C  1 

1  1 

1 

1 

n 

1  .  . 

.  n^  . 

q-1 
.  n^ 

• 

• 

. 

q-1 
.  n^1 

1 

• 

• 

• 

0 

1 

1 

• 

• 

• 

1 

1 

• 

1  .  o 

q-1 

• 

q-1 

.  n 

Proof:  Partition  S  as 


•  ' 
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where  Sq  is  a  (p-q) x (p-q)  upper  triangular  matrix  with 
l's  on  the  diagonal.  Then 


n 

I 

V=1 


Cn  VQ 

b0  blb2 


V 


0 


Divide  the  nonzero  portion  of  Sn  into  four  parts  and 
prove  the  lemma  for  each  part. 


Part  I:  0  <  i  <  j  <  p-q.  Note  that  is  in  the  upper 

right  hand  corner  of  a  ( j-i+1) x  ( j-i+1)  diagonal  block  B 

of  S.  Since  B  is  upper  triangular  with  l's  on  the  diagonal 

Bn=0(n^  1).  Furthermore  (Sn)^j  =  (Bn)^j  w^ich  proves  that 

(Sn)  .  .  =  0(n=>  1)  as  stated  in  the  lemma, 

l  j 

Part  II:  i  >  p-q  and  p-q  £  j  <  p-1.  The  hypotheses 
of  Lemma  2.2  are  satisfied  with  D  =  D.  Therefore  we 
are  entitled  to  use  (2.20),  which  states  that  d^eX^ 
where  X.  =  (S-I)P_1_:IX  ,  and  X  .  is  the  subspace 

j  p— JL  p  — J. 

spanned  by  the  generalized  eigenvectors  of  S  correspon¬ 
ding  to  1.  Thus  there  exist  vectors  x.eX.  such  that 

— _L  JL 


Sn6_j  =  Sn(6_j-dj)  +  Snd  j 


=  Sn  ( 6  .  -d . )  +  d.+  nx.  .  +  ...+ 
~3  ~3  ~3  -D-1 


n 

(J-p+qJ 


x 

-p-q 


+  some  vector  in  X 


p-q— 1  * 


' 


By  part  2  of  the  1st  q-order  condition 


(Sn  ( — j-— j  )  )  i  =  0  for  i  -  P-< 


From  (2.21) 


X  =  <6n  ,  6,  , .  .  .  ,  6  ,  > 

p-q-1  —0—1  — p-q-1 


Hence  for  i  >  p-q 


(Sn6  . )  .  =  0(n3"P+<3) 
— j  1 


Part  III:  i  >  p-q  and  j  >  p-1.  Recall  from  the  proof 

of  Lemma  2.2  that  (z-l)p  is  an  elementary  divisor  so 

that  by  the  q-root  condition  6_.  =  x  +  w  where 

3  P 

x  eX  ,  and  Snw  =  0(nq_1).  Then  there  exist  vectors 
-P-1  P-1  “ 

x.eX.  such  that 

—l  l 


Sn5 .  =  Snw  +  x 


—J 


~  ,  +  nx  „+...+ 

-P-1  -p-2 


n 


Iq-iJ 


X 

~ p~q 


+  some  vector  in  X 


p-q-1 


Therefore  for  i  £  p-q 


(Sn6.)i  =  0(nq  X)  , 


and  so  S^1  is  as  illustrated, 
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Part  IV:  i  <  p-q  and  j  >  p-q.  Consider 

l  Sq^S^"1  =  nO  (Sq  )  0  (S1 )  O  (S2 ) 
v=l 

where  O(Sq)  and  are  as  given  in  the  statement  of 

this  lemma,  and  0(S1)  is  a  (p-q)x(k+q)  matrix  of  l's.  □ 

Corollary: 

| | H—1Sn | |  <  ksh1_p  for  some  constant  kg  depending 
only  on  S.  0 


The  above  corollary  can  be  strengthened  because 


actually 


H-Va.  =  o(h-minli'p"1)>  . 


Therefore  there  exists  a  constant  k^  such  that 


N 


—  I  s  -  £ 


n=0 


— n  ‘  ■  h 


where 


a.  ,h 


P  _ 


i  ,  -1  ,1-p 

col(a0,h  alf...,h  a 


,  1“P  v 

P-1'-"'  ak+p-lJ 


This  bound  suggests  that  R  could  be  measured  by  the  norm 

vi.  ,,P 


. 
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which  is  independent  of  the  method.  Because  this  norm 
is  cruder  than  |.|g,  it  is  unsatisfactory  for  proving 
theorems  in  this  chapter. 

The  function  $  resembles  the  increment  function 
of  a  one-step  method,  and  therefore  it  is  expected  that 
for  a  stable  method,  $(t,a)  should  be  required  to  satisfy 
a  Lipschitz  condition  in  a. 

Lemma  3.2: 

If  a  multivalue  method  satisfies  the  1st  q-order 
condition, then  for  any  Lipschitz  problem,  there  exists 
a  constant  L  such  that 


$(t,a2)  -  S^a1)  |  <L||a2-a1||H 


2  1  2  1 

Proof:  Let  e  =  a  -a  and  e  =  $(t,a  )  -$(t,a  ).  Consider 


three  cases. 


2  1 

Case  I:  The  corrector  only  method,  e  =  <j)  -4>  where 


v=l ,  2 


Because  f  is  uniformly  Lipschitz 


whence 
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(3.2) 


e  |  < 


X  UL.Hh-V 


SH 


i=0 


—  1  H 


p-q  i. 


1  -  l  |  Z.  |L.h 

i=0  ?!  1  1 


p-i 


|  |h_16_TsH|  |  is  bounded  because  .  =  0  for  j  <  i  <  p-q. 
By  the  definition  of  corrector  only  methods,  the  deno¬ 
minator  of  (3.2)  is  no  smaller  than  n .  Therefore  there 
exists  an  L  such  that  |e|  <  L| |e| |^. 


M 

Case  II:  The  P  (EC)  method  with  1  11.  Let  e 

-  P 


(m) 


♦fa)  "  ^  '  and  e 


(m) 


,  x  =  a^  v  -  a^  x  where 
(m)  — (m)  -(m) 


(3.3) 


*  (0) 

Q* 

1 

X 

II 

p!  (AaV 

>p  ' 

u 

♦On)' 

:  (1- 

V  (0) 

+  —  I 

UM  j  =  l 

V 

„  V 

,  ~  hp 

.  V 

-(m) 

=  Sa 

+  Z  — r 
--  p! 

• 

2 

■e- 

m-D. 


v 


(j-D 


Using  the  Lipschitz  condition, 


<  j-D 


(j-1)  1  ‘H 


so  that 


u  m 

(3-4)  e(m)=  (1"S^)£(0)  +  .lx  0(l  I -(j-1)  I  1 H5 


~  hp 

e  ,  %  =  Se  +  Z  — r  £  /  \ 
—  (m)  —  —  pi  (in) 


Also 


• 
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Because  S . .  =  0 
ID 

and  so 


for  j  <  min(i,p-q) , 


is  bounded. 


(3.5) 


=  0(l  ie|  j  H>  +  0(hq"(m))  • 


By  induction  on  m  it  can  be  shown  from  (3.4)  and  (3.5) 
that 


(3.6)  e  =  e  =  0  ( |  |  e  |  |  R)  +  O(hq£^0j) 

Note  that  6TS  =  (1  -  £  )  6T  A  whence 
-p  P  -p 


A  .  =  0  for  j  <  p-q 
PD 


and 


1 I h  p6TAH| | 

u1 

i 

ll 

II  P 

By  (3.3) 


(0) 


p!  h  p5TAe 
-P  - 


=  0  (h“q|  |e|  |R) 

This  result  with  (3.6)  implies  that  £  =  0(||e||H). 

M 

Case  III:  The  P  (EC)  method  with  =  1.  Using  the 
■  P 

notation  of  Case  II,  there  is  considerable  simplifica¬ 
tion  : 


(m) 


=  F  (t ,  a 


— (m-1) 


)  -  F(t,a 


-(m-1) 


) 
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G  /  -»  \  Se  ^  *  G  /  **v  • 

—  (m-1 )  —  —  p!  (m-1) 


Using  the  Lipschitz  condition. 


(m) 


p-q  .  , 

<  T  L.  ^-1  (Se)  .  +  Z. 
i=0  1  h1  1  1 


h£ 

p!  £ (m-1) 


Because  Z^  =  0  for  p-c  <  i  <  p-q,  we  have 
e  (m)  =  °U  ls£l  lH)  +  0  (h  e  (m-i)  ) 

(3.7)  =  0  ( |  |  e  |  |H)  +  0(hce(m_1))  . 

By  the  predictor  condition  for  r  =  1, 

A  .  =  0  for  j  <  p-Mc 
PD 

whence 

|  |  h_P6_pAH  |  |  =  0  (h-Mc) 

By  (3.3) 

E  (0)  =  p!  h-P^ 

(3.8)  =  0(h"Mc| |e| |H)  . 


Together  (3.7)  and  (3.8)  imply  e  =  e =  0 ( | | e | | R) . 


*  '  i 


(€.£)  yj.€ 
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The  following  lemma,  proved  in  Miller  (1968:21), 
is  the  discrete  counterpart  to  Gronwall's  inequality. 

It  is  useful  for  solving  inequalities  which  arise  in 
stability  proofs. 


Lemma  3.3: 

Let 

n-1 

w  <  y  +  T  v.w.  for  n=0,l,2,... 
n  j=0  3  => 


where  y,  w  ,  and  vn  are  nonnegative.  Then 


y  exp 


n-l 


l 

j=0 


□ 


Theorem  3.2: 

A  multivalue  method  which  satisfies  the  q-root 
condition  and  the  1st  q-order  condition  is  q-stable. 


Proof :  Consider  a  Lipschitz  problem,  and  set  e_n  = 

a  -aX.  Then 
— n  — n 


e  =  Se  t  +  Z  — r  e  +  r 

— n  —n-l  —  p!  n  — n 


where  en  = 

exists  a  constant  L  such  that 


and  by  Lemma  3 . 2  there 


(3.9) 


H 


* 


V 


Solving  the  difference  equation  gives 


„n  ,  r  on_ j  /n 

en  =  s  +  I  S  U  jjr  e.  +  r. 


Therefore 


,  p  n  .  n 

e  |  |  w  l  I £  •  I  I  |sn  D^|  L  +  II  £  S*  3  r  .  ,  ,  „ 

— n  1  1  H  ”  p !  .  ^ ,  1  3  1  11  — 1  1  H  -A  — 3  1  1  H 

^  3=1  J  3=0  J 


From  the  corollary  to  Lemma  3.1  it  follows  that 


Sn  I H  *  kshl  P!  11 


This  inequality  and  (3.9)  imply  that  there  is  a  constant 
k^  such  that 


n 


sJIh  s  hki  .I.llaj-i11-  +  IR 

3=1 


3-±'  'H  '  1  — 1 S  * 


Apply  Lemma  3.3  with  w  =  | | e n | |H,  y  =  |r| g,  and  vn  = 
hk^  to  get 

||e^||H  <  |  R  |  s  expCnhk^j^)  <  |R|g  exptk^  . 


Together  with  the  corollary  to  Lemma  3.1,  this  proves 
q-stability.  □ 


Section  3.3  Proving  rth  Order  q-consistency 

The  proof  of  rth  order  q-consistency  assumes  the 

g 

rth  q-order  condition,  and  so  the  choice  of  a  (t) =  DY(t) 


I  !  S  II  T 

automatically  satisfies  |  |a  (t)-a(t)  |  | H  =  0  (h  ) .  To 
show  that  I  I Sn  (aS  (0 ) -a (0 )  I  I u  =  O  (hr)  ,  we  first  prove 
the  following  lemma,  which  has  the  q-root  condition  as 
an  additional  hypothesis. 


Lemma  3.4: 

If  a  multivalue  method  satisfies  the  q-root  con¬ 
dition  and  the  rth  q-order  condition,  then 


(Sn(D-E))i  . 


0(nj"r‘min(i'p_q))  . 


Proof:  Let  u.  = 

-  —j 

of  values  for  j . 


(D-E)e.,  and  consider  four  sets 


Case  I :  j  <  r.  Then  u..  =  0_. 

Case  II:  r  <  j  <  r+p-q.  (u . ) . =  0  for  i  >  j-r.  From 

3 

Lemma  3.1  it  can  be  seen  that 


(Snu  . )  .  =  0  (n-,-r_1) 
~3  i 


Case  III:  r+p-q  <  j  <  r+p-1.  By  the  rth  q-order  condi¬ 
tion 

(Snu . ) .  =  0(n^~r~P+q)  for  i  >  p-q  . 

— j  l  - 


Partition  S  as  in  Lemma  3.1  and  likewise  partition 
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u .  = 


u 


u 


Thus 


n 


Snu. 

~3 


s^u1  +  i  sj-vrv 1 

v=l 


on  2 

S2u 


where  S^u2  =  O (n^  r  p+g) .  The  size  of  Sq  is  known  from 
Lemma  3.1,  and  so  it  is  easy  to  show  that 


(Snu  ).  =  0(n3-r-min(i'P-<3))  . 


Case  IV:  j  =  r+p-1.  From  Lemma  3.1,  it  is  obvious  that 


Snu  ,  ,  =col(0(n^  ■*■),...  ,0  (n^  ^ )  ,  .  .  .  ,0  (n^  ^)  )  . 

— r+p-l 


□ 


Theorem  3.3: 

If  a  multivalue  method  satisfies  the  q-root  con¬ 
dition  and  the  rth  q-order  condition,  then  it  is  q- 
consistent  of  order  r. 


Proof:  With  as  (t)  =  DY (t )  the  first  requirement  of  q- 

consistency  follows  from  part  1  of  the  rth  q-order  con¬ 
dition.  Applying  Lemma  3.4  with  h  1  substituted  for  n 
gives 


■ 
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[Sn  (aS  (0 )  -a  (0 ) )  ]  .  =  (Sn(D-E)Y(0))  . 
—  —  1  —  1 


r+p-1 

l  (Sn (D-E) ) . .0 (hj ) 
j  =  0 


=  O  (hr+rain  ) 


whence 


|  |sn(ab  (0)-a(0)  )  |  |R  =  0(hr)  , 


thus  proving  the  second  part  of  q-consistency .  Recall 
that 


(3.10)  PY(t-h)  =  Y(t)  +  0(hr+P)  . 


Use  this  and  Lemma  2.1  to  deduce  that 


Sa  (t-h)  =  SDY(t-h) 


=  DY  (t)  -  l  eTY  (t)  +  0(hr+P) 

-  - p-  - 


(3.11) 


=  ab  (t)  -  £ 


hP  (p) 


.  y 


(t)  +  0(hr+P) 


By  definition 


(3.12)  as(t)  =  Sas  (t-h)  +  £  ^  $(t,aS(t-h))  , 


and  so  we  need  only  show  that 


$ (t ,as  (t-h)  )  =y(p)(t)  +  0(hr)  . 


■ 


' 
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Case  I:  The  corrector  only  method.  We  have  $ (t ,aS (t-h) ) = 
<p  where  <f>  =  F(t,aS(t)).  Also  y  ^  (t)  =  F(t,a(t)),  and 
by  using  the  Lipschitz  condition, 

(3.13)  U-y(p)(t)|  s  Pyq  L.  if  la? (t)-a.  (t)  I  . 

i=0  1  h  1 


From  (3.11)  and  (3.12) 

as(t)=as(t)+i^  (*-y(p)(t))  +  0(hr+P)  . 

P  • 

For  i  <  p-q 


af(t)-a.  (t)  =  £.  ~  (<j>-y (p)  (t)  )  +  0  (h1+r)  . 

JL  -L  JL  |J  • 


Substitute  this  into  (3.13): 


(p-y (p)  (t)  I  <  ~r  L  i  ±4-  \z.  (^-y^'  (t)).+  o(hA'ri)  | 

i=0  h1  p* 


p-q 


i! 


i+r, 


which  gives 


p-q 

l  Li 

i=0 


\l±\)  |  4>-y  (p)  (t) 


=  0  (hr )  . 


Since  the  first  factor  is  no  smaller  than  n  #- 

<j)  =  y  ^  (t)  +  o  (hr) 

M 

Case  II:  The  P (EC)  method  with  ^  1 .  For  this  case 

p 

$ (t,as (t-h) )  =  where 

(3.14)  <j>(0)  =  h“p  p!  (Aas  (t-h)  )  , 
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u  i  m 

mx  ,  .1  r  /  t  vin-]. 


(3-15)  ♦(«,)-  (1-^>  *(0)  +  ^  .1  (1-V  Jp(t'*(j-1))' 


a(m)  =  %S  <t-h)  +  I  £r  * 


p!  T  (ra)  * 


Let  £(m)  =  <!>  <m)  “  y(P)(t)  and  =  a(n,l  -  a(t)‘  Note 


—  (m)  ^(m) 


that 


(  \  ^  /  \  t  in  • 

(3.16)  y(p)  (t)  =  (1-  -2)y  (p)  (t)  +  -f-  \  (l-£  )m_:iF  (t,a  (t)  ) 

M  M  j=l  p 


Subtract  (3.16)  from  (3.15),  and  use  the  Lipschitz  pro¬ 
perty  to  get 


u 


(3-17)  £  (m)  =  (1  -  5^  £  (0 )  +  1,  0(l  l-(j-l)  l  'H‘ 


m 

l 

j=i 


Also  by  (3.11) 


— (m)  =  S^(t"h)  +  I  |r  ^  Cm)  ■  2-(t) 

,  P 

=  aS(t)  +  0(hr  P)  +  l  “y  £  (m)  -  a(t) 

“  i  S  £  (m,  +  H° 

and  so 

(3.18)  I  l2.(m)  I  lH  =  0(hr)  +  0(hq£(m))  • 

Write  Lemma  2.1  as 

(I-I  ST)AD  =  (I-£  6T)EP  +  (D-E)P 

- p  - p 


' 
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whence 


(3.19)  6TAD  =  eTP  +  (l-£  )  1  ST(D-E)P 
~P  “P  P  ~P 


so  that  using  (3.10) 


6TADY(t-h)  =  Y  (t)  +  0  (hr+P  q)  . 

“P  -  P 

It  follows  from  (3.14)  that 

e(0)  =  h"p  p!  ^ADY(t-h)  -  y(p)(t) 

=  0(hr-3)  . 

This  result  with  (3.17)  and  (3.18)  implies  that 


u  _ 

(3.20)  e(m)  =  (1-  ^)e(0)  +  0(hr)  =  0(hr  q)  , 

and  consequently  <J)  -y^  (t)  =  e  =  0(hr). 


(M) 


M 

Case  III:  The  P  (EC)  method  with  Z  =1.  The  expression 

P 

for  e  is  very  simple: 


e(m)  =  F(t'2-(m-l)>  '  F(t'*(t))  ' 


and  as  in  case  II 


r  hF 

e ,  =  HO(h  )  +  Jl  n-  e,  \ 

—  (m~ 1 )  —  -  p!  (m-1) 


The  Lipschitz  property  enables  us  to  write 
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p-q 


|  e  ,  J  <“  l*  L.  ii-lo  (h1+r)  +  £.  ~  e  ,  ,  . 

i=o  1  h1  1  p  (m-1) 


and  because  =  0  for  p-c  <  i  1  p-q  , 


e  (m)  “  0(hr>  +  0(hCE(n.-l))  * 


Hence 


(3.21)  c(m)  =  0(hr)  +  O(hmce(0))  . 


Part  3  of  the  rth  q-order  condition  states  that 


(AE)  .  =  A  .  =  P  .  for  j  <  r+p-Mc  . 
PD  PD  PD 


Thus  A  .  =  0  for  i  <  p-Mc.  Let  j  <  r+p-Mc  so  that 
Pi 

(D-E)^^  =  0  for  i  >  j-r.  Since  j-r  <  p-Mc, 


(A (D-E ) )  .  =  0 

J 


for  j  <  r+p-Mc  , 


and  hence 


(AD)  .  =  P  . 
PD  PD 


for  j  <  r+p-Mc  , 


which  implies 


(ADY(t-h))  =  Y  (t)  +  0(hr+p  Mc) 


By  (3.14) 


. 
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e  =  h  Pp!  (ADY(t-h))  -  yv*"(t) 


(P) 


(0) 


„  /v. r-Mc, 
=  O  (h  ) 


This  result  with  (3.21)  implies  that  e 


(M) 


0(hr).  □ 


CHAPTER  4 


ASYMPTOTIC  BEHAVIOUR  OF  THE  ERROR 


In  this  chapter  asymptotic  estimates  of  the 

global  error  are  found  for  corrector  only  methods 
M 

and  P  (EC)  methods.  It  is  shown  that  the  global  error 

e  defined  by  e  =  a  -a (t  )  satisfies  an  equation  of 
— n  — n  — n  —  n 

the  form 

H_1e  =  hre (t  )  +  0(hr+1) 

— n  —  n  — 

provided  that  f  (t ,y , . . . ,y )  is  sufficiently  smooth 
and  there  is  no  round-off  error.  Here  the  function 
e  (t)  is  the  magnified  error  function,  which  depends  on 
t  alone. 

The  principal  results  of  this  chapter  are  con¬ 
tained  in  Theorem  4.1.  Because  the  statement  of  this 
theorem  is  long  and  the  proof  is  tedious,  definitions 
of  various  quantities  appearing  in  the  theorem  are 
presented  here.  The  definitions  of  the  matrices  D  and 
E  are  changed  for  this  chapter  only.  Instead  of  D  and 
E  being  (k+p)x(r+p),  they  are  (k+p) x  (r+p+1) ,  and 
similarly,  has  r+p+1  components  for  this  chapter. 
Furthermore,  the  Pascal  triangle  matrix  P  is  redefined 
to  be  a  (r+p+1) x  (r+p+1)  matrix.  Lemma  2.1  clearly  holds 
for  these  enlarged  matrices: 
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. 
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(4.1)  SD  =  DP  - 


T 

*  %>  p 


Also  define 


Y(t)  =  col(y(t)  ,hy'  (t)  ,  .  . .  ,hr+py  (r+p)  (t)/(r+p)  !) 


and 


3f 


, . . . ,y  (p  q)  (t) ) ,  i=0,l, . . . ,p-q  , 


=  — nr 

dy 

/  >  \ 

assuming  3f/3y  1  exists. 

The  following  lemma  is  used  in  the  proof  of 
Theorem  4.1. 


Lemma  4.1: 


Let  3f/3y^  have  continuous  partial  derivatives 


for  te[0,l],  and  let  32f/3y  (l) 3y  ( j }  be 
te  [0 , 1] .  Then 


bounded  for 


(4.2) 


F(t,a)  =  y^p)  (t)  +  l  J.  (t)  (h  xi !  a .  -y  ^  (t) )  1 
“  i=0 


p-q 


l  0  (h  1i  !  a  .  -y  '"L;  (t)  ) 
i=0  1 


(i) 


Proof:  By  two  applications  of  the  mean  value  theorem 

p-q  3F (t ,a (t) ) 

F(t,a)  =  F (t ,a  (t)  )  +  l  — 

i=0 


8  ci . 
l 


^ai~ai  ^  ^ 


p-q  p-q  32F(t,a(t)  +0i[a-a(t)]) 


+  i  .K 

i=0  j  —  0 


3  a  .  3  a  . 
i  3 


(ai~ai  (t)  )  (a^a^  (t) ) 


where  0^e(O,l)  depends  on  t  and  a.  The  lemma  follows 

by  noting  that  3/3a.  =  h* 1!!  — ^-r-y  .  □ 

1  8yUJ 

The  following  theorem  gives  the  asymptotic  beha¬ 
viour  of  the  numerical  solution  for  a  special  class  of 
starting  values. 

Theorem  4.1: 

Consider  an  rth  order  q-convergent  multivalue 

M 

method.  For  a  P  (EC)  method  it  is  required  that  if 

i  =1  and  0  <  r+p-Mc  <  p-c  then  A  ,  ,  _  =  0  .  Let 
y (t) sCr+P+^ [ 0 , 1]  be  the  solution  of  a  Lipschitz  problem 
for  which  f  satisfies  the  hypotheses  of  Lemma  4.1.  Let 
the  starting  value 


aA  =  DY ( 0 )  + 


where 


(  T3  ""  1  ) 

6^  '  are  freely  chosen.  Then 


where 


i  p;1 

a‘L(t)  =  DY  (t)  +  l  d. 

i=0 


The  function  6 (t)  solves  the  initial  value  problem 


-V 
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6  (l)  (0)  =  6(jl) 


for  i=0 , 1 , . . . ,p-l  , 


(P) 


p-q  (is 

(t)  =  l  J.  (t)  (6  u;  (t)  + 
i=0 


i! 


(i+r)  !  i  ,  i+r 


D,  ,  .  _y  (l+r)  (t)  )  +  x(t) 


where  x(t)  is  specified  for  each  method  as  follows 


1.  for  a  corrector  only  method  and  for  a  P (EC) 


M 


method  with  £  =  1  and  r+p-Mc  <  p-c 


x(t)  =  0  . 


2. 


M 

for  a  P (EC)  method  with  £  ^  1 

P 


-  *  ,  £  M(l-£  )M  1 

x(t)  =  ■— (i - P - )  (D-E)  x 

(r+p-q)!  £p  uM  p,r+p-q 


J  „  (t)  y  (r+^>*"<^)  (t)  . 


p-q 


3. 


M 

for  a  P  (EC)  method  with  £  =1  and  r+p-Mc >  p-c 

P 


x  (t)  = 


£ 


(p-c) ! 
p-c  p! 


J  (t) 
p-c  ^ 


M 


P 


v,  (AD-EP)  .  ..  x 

(r+p-Mc) i  p, r+p-Mc 


(r+p-Mc)  {t)  _ 


Proof:  Define 


A (t)  =  col (5 (t) ,h6 '  (t) , . . . 5  (p)  (t) ,0, . . . ,0) 

hr  • 


and  aS(t)  =  D  (Y  (t)  +  hrA  (t) ) .  Since  a1(t)=  as (t ) +0 (hr+p) , 


Ii  S  |  | 

a  -a  (t  )  „ .  By  the 

1  — n  —  n  1  1  H  J 

triangle  inequality 


By  q-stability  a  bound  on  the  first  term  is  given  by 
i  I  a  -aX I  I „  .  kn I R I _ 

A  bound  on  the  second  term  is  found  by  using  q-stability 
and  the  fact  that  as  (t-Q )  , as  (t^)  ,  .  .  .  , as  (t^)  is  the  compu¬ 
ted  solution  when 

Rs  =  (as  (tQ)  -a1  (tQ)  ,  as  (t1)  -as  {t±)  ,  .  .  .  ,as  (tN)  -as  (tN)  ) 
where 

,  p 

(4.3)  aS (t) =  SaS  (t-h) +  I  $  (t,aS  (t-h) )  . 

Thus 

1  Ias  (tn)-a^l  lH  ^  k0kshl~PllaS<to)-a1(to)M 

+  k0kShl”P  ^  I  |as(t;.)-as(t.)  |  |  . 
j=l 

Therefore  we  need  only  show  that 

Ss(t)  =  as(t)  +  0(hr+p+1)  . 


••••  *  **  ;c  >  >" l3t  ^  Jj 
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Since  6  (t) £CP+^  [0 ,1]  , 


(4.4) 

P (Y (t-h)+hrA  (t-h) )  =  Y  (t)+hrA  (t) +  0 (hr+p+1)  . 

Use  this 

and  (4.1)  to  deduce  that 

(4.5) 

Sas  (t-h)  =  aS  (t)  -  z  ^  (y(p>  (t)+hr6(p>  (t)  )+0  (hr+p+1) 

P  * 

After  substituting  (4.5)  into  (4.3),  it  becomes  clear 
that  we  need  only  show  that 


(4.6) 

$  (t ,as (t-h) )  =  y (p)  (t) +  hr6 (p)  (t) +  0  (hr+1)  . 

Case  I : 

The  corrector  only  method.  In  this  case 

<Ht,as(t-h))  =  <P  =  F  (t ,aS  (t)  )  . 

From  (4.3)  and  (4.5)  it  follows  that 

Is  (t)  =  aS  (t)  +  l  ^  (4-y  (p)  (t)  )  +  O  (hr+p)  . 

P  • 

In  the  proof  of  Theorem  3.3  it  was  shown  that 


and  thus 

4>  =  y  (p)  (t)  +  0  (hr) 

for  i  <  p-q 

h_1i la?  (t)  =  y(l)(t)  +  hrw± (t)  +  0(hr+1) 

where 


where 


w.  (t)  =  6  (l>  (t)  +  -(-■  D.  .  y <i+r)  (t)  . 

i  (i+r) !  1,1+^ 


Apply  (4.2)  with  a  =  a  (t)  to  get 


4)  =  y  (P)  (t)  +  "  l  "  J  .  (t)  (hrw.  (t)  +  0  ( h r'r-L) )  +  O  (hzr) 

i=0  1  1 


p-q 


r+l. 


2r, 


=  y(p)  (t)  +  hr6(p)  (t)  +  0(hr+1) 


thus  showing  (4.6). 


M 

Case  II:  The  P  (EC)  method  with  £  ^  1 .  In  this  case 

$  (t,as  (t-h)  )  =  4>  where 


=  h  pp!  (AaS  (t-h)  )  , 


I  d“An)m  Dp(t,ar.  n) 


u  ,  m 

^(m)  “  (1  "  u^)C|>(0)  +  j£x  *  ^'-(j-l) 


-(m) 


s  ~  hp 

=  Sas  (t-h)  +  Z  cj) 


p!  (m) 


Let  e(m)  =  *(m)  '  y(P>  (t)<  By  (4-5)' 


*(m)  =  ^{t)  +  i  h  (hm)“y(P>  (t))  +  0<hr+P) 


=  aS  (t)  +  l  |j-  e  (m)  +  O  (hr+P) 


so  that  for  i  £  p-q 
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(4.7)  h*~ii I  (2L(m)  )j_-  y  (i)  (t)+hrwi  (t)+hp  1  ■  .-1-  e 


.  i !  &  . 


(m) 


+  0(hr+1) 


From  (3.20) 


u 

(4 . 8  J  e(m)  =  U-^o)  +  O  (hr)  =  0(hr  q)  . 


Apply  (4.2)  with  a  =  a^  to  get 


p-q 


r(t,a(m))  =  y'p  <t)  +  hr  I  ji(t)wi(t) 


+  J  (t)  (p  q)-!  hqI  e,  .  +  O  (hr+1) 
p-q  p!  p-q  (m) 


Use  this  result  to  obtain 


M 


(M)  u. 


=  TT-  l  (l-^)M_:lF(tra{j_i)) 


M  j=l 


(4.9) 


(„\  r  p-q 

=  y  P  (t)  +  hr  l  J. (t)w. (t) 

i=0 


l  M 


+  iE^hq  -Eia  .I.(l-^)M-jc(j.1)+0(hr+1) 


p-q 


UM  j=l  P 


Recall  (3.19) 


STAD  =  cTP  +  (l-£  )  1  6^(D-E)P 

-P  -P  P  “P 


Using  (4.4) 


(0) 


h  Pp !  (AaS  (t-h)  )  -  y  (p)  (t) 


=  In  PnUT 


(P) 


h  ppl  6  AD  (Y (t-h)  +  h  A (t-h)  )  -  y  vjy'  (t) 


-P 


h-pp!  (et-fl-I  )  (D-E)  )  (Y(t)+hrA(t))-y(p)  (t) +0  (hr+P+1) 

P  P  P 


=  h‘ 


p!  (D-E)  ,  0  « 

‘-q  _ _ P  r  r+p-q  (r+p-q) 


(r+p-q) 1  (1 -l  )  Y 

hr 


(t)  +  0  (hr+q+1)  . 


Substitute  this  into  (4.8),  and  then  put  the  result  into 
(4.9).  After  considerable  algebraic  manipulation,  it 
turns  out  that 

♦  (M)  =  y <p>  (t)  +  hr«(p)(t)  +  0(hr+1)  . 

M 

Case  III:  The  P  (EC)  method  with  %  =1  and  r+p-Mc >  p-c. 

p 

In  the  proof  of  Theorem  3.3,  (3.22)  clearly  holds  for  the 

augmented  D : 

(AD)  .  =  P  .  for  j  <  r+p-Mc 
P3  P3 

whence 

(6TAD) .  =  ((eT  +  aeT,  M  )P) .  for  j  <  r+p-Mc 
— p  j  — p  —r+p-Mc  j 


where 


a  =  (AD-EP) 


p, r+p-Mc 


Therefore 


6  ADY (t-h) 

-p  - 


y(P)  (t)  +  (r+p-Mc)  (t)+0(h 

p !  2  (r+p-Mc) !  u 


r+p-Mc+1 


■ 


and  from  (4.4) 


e(0)  =  h  pp!  (ADY(t-h)  )  -  y(p)(t) 


(4.10) 


=  hr'Mc  ,  x  y(r+P’Mc)  (t)+0(hr’Mc+1) 

(r+p-Mc) I  J 


From  (3.21),  we  know  that  for  m  <  M 


e(m)  =  O(hr)  +  O(hmcE{0))  =  0(h 


r+  (m-M)  c 


) 


Thus  for  i  <  p-q 


.  i  I  it 


pi  (m-1) 


0(hr+(m-M)o+l) 
,  c  <P-c)!Vc 


(m-1) 


0 


if  i  <  p-c 
if  i  =  p-c 

if  i  >  p-c 


It  follows  from  (4.7)  with  1  <  m  <  M  that 


h  1i I (a (m_1) ) ±  =  y (l) (t)  +  hrw± (t)  +  O (hr+ (m  M)c+1) 


c  (p-c)  I 

pi  A'p-cc  (m-1) 


h~  ^  T#  -  i  n  if  i  =  p-c 


0 


otherwise 


Apply  (4.2)  with  a  -  a-.(m-i) 


to  get 
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(m) 


=  F<t'2.(m_i))  -  y(p)  (t) 


hr  Pyq  J.  (t)w.  (t)  +  hcJ  (t)  (p  e  ,  . 

i  i  p-c  p!  p-c  (m-1, 


+  0(hr+(m-M)c+1)  +  0[h2<r+(m'M)c)]  . 


(4.11)  =  hr  £  J.  (t)w.  (t)  +  hcJ 

i=0  p  p* 


p !  ^p-c£ (m-1) 


+  0(hr+(m"M)c+1) 


because  r+(m-M)c  >  r+(l-M)c  >  0.  Notice  that  if  m  <  M 
then  the  first  term  of  (4.11)  can  be  absorbed  into 
0(hr+(m-M)c+lK  Therefore  using  (4.10)  and  (4.11)  to 

solve  for  £ (M)  ?ives 

e  (M)  =  hr6(p)  (t)  +  0(hr+1)  , 


which  proves  (4.6). 


M 

Case  IV:  The  P(EC)  method  with  =  1  and  r+p-Mc  <  p-c 
“  "  t' 

Because  of  the  restriction  on  the  method 


A  .  =  0 
PJ 


for  j  £  r+p-Mc  , 


which  implies 


e  =  h  pp! (Aas (t-h) )p  -  y  ^  (t) 


=  0(hr-Mc+1)  , 


y  —  C+  1  . 

and  so  =  0  (h  ).  Thus  for  l  £  p-q 


It  follows  from  (4.7)  that 


h'h!  (a(M_1))i  =  y* * * * * * * (i) * * * * * *(t)  +  hrWi(t)  +  0(hr+1) 

from  which  it  is  straightforward  to  get 
e,„,  =  hr6(p)(t)  +  0(hr+1)  . 


Corollary: 

Assume  that  the  method  is  strongly  stable  and 

that  f  =  f(t).  Let  0  <  e  <  1  be  given.  Then  for  a^ 

a  (0 )  and  r_±  =  ...  =  ^  =  0,  we  have 


I  d.  T 

a  =  a  (t  ) +  HO (h  )  uniformly  for  t  e[e,l] 

— n  —  n  —  n 


where 


(i)  _ 

0 


Proof:  For  the  differential  equation  y  P^  =  f (t)  it 

is  easily  seen  that 

a  =  aX  +  Sn (a (0)  -  a1(0))  . 

— n  — n  —  — 

1  v-4- 1 

By  Theorem  4.1,  a*  =  a  (tR)  +  HO(h  )  ,  and  hence  we 
need  only  show  that 
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(4.12)  Sn  (a  (0)  -  a1  (0) )  = HO  (hr+1)  for  n  >  eN. 

To  verify  (4.12),  it  is  sufficient  to  show  for  i  =0,1,..., 
p-1 

(4.13)  Sn[  (D-E)ci+r-  cL  n_i(Sn(D-E))0|i+r]=  HO(h1_1) 

for  n  ^  eN.  Let  be  defined  as  in  the  proof  of  Lemma 
2.2,  and  let  W  be  the  space  spanned  by  all  generalized 
eigenvectors  corresponding  to  eigenvalues  other  than  1. 
Because  the  method  is  strongly  stable,  dim  =  j+1. 

__  _ I 

Furthermore,  if  weW  then  S  w  =  0(N  )  for  n  >  eN .  By 

the  rth  q-order  condition,  if  i  <  p-q  then 

(E>-E)£i+r  £<60,61,...,6i>  =  xi  =  <d0,d1,...,di>  . 

The  coefficient  of  ch  in  (4.13)  has  been  chosen  so  that 
the  expression  in  square  brackets  is  a  vector  in 
With  this  information  it  is  straightforward  to  verify 

(4.13) .  Using  Lemma  3.4,  it  is  not  difficult  to  show 

that  if  i  >  p-q  then  (D-E)  £i+reXi +  W.  The  coefficient 
of  d^  has  been  chosen  so  that  the  expression  in  square 
brackets  is  in  Xi_1+  W.  Once  again,  the  verification 
of  (4.13)  is  straightforward.  □ 

Results  of  Henrici  (1962:255)  cause  us  to  believe  that 
this  corollary  is  also  valid  for  f  =  f (t ,y , . . . ,y ) . 


. 


CHAPTER  5 


VARIABLE  STEPSIZE 


Here  let  us  consider  the  general  grid 

0  -  t0  <  ti  <  ...  <  fcN  -  1. 

To  measure  the  fineness  of  the  grid,  define 


h  = 


max 

1  <  n  <  N 


h 


n 


where  h  =  t  -  t  , .  For  this  chapter  the  dependence 
n  n  n-1 

of  various  quantities  on  the  stepsize  h^  is  made  expli 
cit  in  the  notation. 

Let  a  =  h  /h  ,  where  hn  =  h, ,  and  define 
n  n'  n-1  0  1 

C  =  diag(l,a  , . . . ,ak+p_1) .  Notice  that  we  have 
n  n  n 

a(tn_i;hn)  =  cn-(tn-l;hn-l)  *  Th±S  ec3uality  is  the 
rationale  for  the  interpolatory  technique  of  varying 

stepsize . 


Definition  of  variable  stepsize  multivalue  method: 
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How  the  grid  is  chosen  does  not  concern  as  long 
as  it  satisfies  certain  restrictions.  For  this  chapter 
it  is  assumed  that  there  exists  a  positive  constant  A 
independent  of  h  such  that 

A  <  6  <1 

n 

where  the  relative  stepsize  0n  =  hn/h.  In  addition,  for 
Section  5.1  it  is  required  that  there  exists  a  positive 
constant  V  independent  of  h  such  that 
N 

y  I  0  -  0  ,  I  £  V  , 

L.  1  n  n-1 1 
n=l 

and  for  Section  5.2,  which  discusses  (k+1) -value  Adams 
methods,  it  is  required  that  changes  in  stepsizes  do 
not  occur  more  often  than  every  k-1  steps.  These  res¬ 
trictions  on  the  stepsizes  are  weaker  than  those  con¬ 
sidered  by  Tu  (1972)  for  the  interpolatory  technique 
of  changing  stepsize. 

In  Section  5.1  it  is  shown  that  the  1-root,  con¬ 
dition  and  the  rth  1-order  condition  are  sufficient  for 
variable  stepsize  convergence  of  order  r,  and  in 
Section  5.2  (k+1) -value  Adams  methods  are  shown  to  be 

kth  order  convergent. 

Section  5,1'  Multivalue  Methods 

For  this  section  it  is  required  that  there  exists 
a  positive  constant  V  independent  of  h  such  that 


■ 
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N 


l  I  6 

n=l 


n 


e 


n-l 


<  V  . 


When  the  stepsize  is  varied  subject  to  this  restriction, 
q-convergence  is  not  preserved  unless  q  =  1.  There¬ 
fore  only  the  case  q  =  1  is  considered  here. 

In  order  to  show  the  significance  of  the  restric¬ 
tions  on  the  stepsize,  let  us  suppose  that  the  stepsize 
is  selected  by  means  of  a  function  0 (t) : 


t  ,  +  h0  (t  , ) 
n-l  n-l 


(Ignore  the  minor  problem  of  getting  tN=  1.)  Then  it 
is  required  that  0 (t)  be  bounded  away  from  zero  and 
that  0 (t)  be  of  bounded  variation. 

The  meaning  of  variable  stepsize  is  sufficiently 
general  to  include  practical  schemes  for  choosing  step- 
size.  Certainly  any  worthwhile  scheme  should  satisfy 
the  two  requirements  involving  A  and  V.  A  more  restric¬ 
tive  kind  of  stepsize  selection  is  considered  by  Tu 
(1972),  who  requires  that  1 0  -  0  ,  |  <  h  0''.  He  shows 

that  if  h  is  controlled  so  that  the  local  error  estimate 
n 

is  equal  to  e  or  ehn  then  this  requirement  is  met.  How¬ 
ever,  for  some  practical  schemes,  0R  does  not  satisfy 
IB  -  6  ,  I  <  h  0'  for  some  constant  G''.  For  example, 

the  computer  program  DIFSUB  (Gear  (1971) )  does  not  per¬ 
mit  increases  of  less  than  ten  percent  in  the  stepsize, 
and  so  | 0n  -  en_1 I  ^  A/10  whenever  0R  >  6n_1. 


* 
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The  theory  of  previous  chapters  extends  to  the 
variable  stepsize  case.  The  definitions  and  theorems 
are  restated  for  variable  stepsize  in  a  slightly 
abbreviated  form.  Let  us  adopt  the  abbreviated  nota- 

TD""  1 

tion  || . | |  =  | | . || H  where  Hn=  diag (1 ,hR, . . . ,hn  - - - 

hn_1)- 

Definition  of  rth  order  convergence: 

A  multivalue  method  is  convergent  of  order  r  if 
for  any  Lipschitz  problem  with  solution  y(t)eC  ^[0,1], 


I  I  n  =  0  (h  ) 


whenever  r^  is  a  function  of  h^  satisfying 


and  r ^  =  r 2  =  ...  =  ^  -  0_* 


□ 


With  q  =  1,  it  is  convenient  to  use  definitions 
of  stability  and  consistency  which  are  more  conventional 
than  the  definitions  of  1-stability  and  1-consistency 
given  in  Chapter  3 . 


Definition  of  stability: 

A  multivalue  method  is  stable  if  for  any  Lipschitz 
problem,  there  exists  a  constant  kQ  such  that 


. 
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□ 


In  order  to  show  stability  for  variable  stepsize,  it  is 

"“In 

necessary  to  have  | | H  S  h| |  bounded,  whereas  for 
fixed  stepsize  it  is  sufficient  to  require  only  that 
h^  |  |  H  1Sn  |  |  be  bounded. 


Definition  of  rth  order  consistency; 

A  multivalue  method  is  consistent  of  order  r  if 
there  exists  a  matrix  D  such  that  for  any  Lipschitz 
problem  with  solution  y  (t) eCr+^ [ 0 , 1] ,  as  (t;h)  =  DY (t ;h) 
satisfies 

1.  ||aS(t;h)  -a(t;h)||H  =  0(hr)  . 

2.  aS(t;h)  =  aS(t;h)  +  0(hr+^)  .  □ 

Theorem  5.1: 

A  multivalue  method  which  satisfies  the  1-root 
condition  and  the  1st  1-order  condition  is  stable. 


Proof :  Consider  a  Lipschitz  problem,  and  set  e^=  a^-a^ 


e 

— n 


SC  e  , 
n— n-1 


r 

— n 


Then 
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x 


where  en  =  *  -  4>  ,  and  by 

Lemma  3.2  there  exists  a  constant  L  such  that 


e  <  L  C  e  , 
n  1  1  11  n— n-1 1  1 n 


Write 


(5.1)  e  =Se  ,  +  g  +r 

v  '  — n  -n-1  — n  -n 


where 


2n 


h 

S  (C  -  I)e  ,  +  £  -?■ 
n  — n-1  —pi  n 


Solving  the  difference  equation  (5.1)  gives 


-n 


l  Sn"jaj  +  l  Sn-^r 

j=l  3  j=0  3 


Hence 


n 


<5'2>  I  I  — n  I  I  n  S  HI  Vs  HjH  H2j"3 


n 


+  y  |  I H“1Sn“^H .  [ I  I  I r . 
j£0 11  n  3 1 1  “3D 


By  using  the  bound  on  e  ,  it  follows  that 


q  II  <  (  I  I H  1SH  -| 

•in  1  1  n  ~  11  n  n-1 


C  -I 
n 


+  -2.1  l^hP”1!  |  L  |  |  H  "'C  H  , 
n ! 1  1 —  n  n  1  1  n  n  n— 1 


,-l. 


— n-1 1 1 n-1 
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Consider  the  (i,i)th  component  of  Cn~  I 


<ciTI,iil  =  l(Ven-l)1  -  1 


< 


A  1  -  1 


-  1  -  A 


l9n  '  9„-l 


Since  | I H  1SH 


—  ,  I  I  ,  I  I £hP  1 1  I  ,  and  I  |H  "C  H  , 

n  n-1 1 1 '  1 n  1 1 n  1 1  n  n  n-1 


-1. 


are  bounded,  there  exist  constants  Cq  and  c^  such 
that 


e  t 

—n-1 


n-1  * 


From  Lemma  3.1  we  see  that  there  is  some  k^  such  that 
|  |  Hn"*"Sn  |  |  <  k^.  Apply  Lemma  3.3  to  (5.2)  with 


v 

n 


-1 


=  ki(co 


e 

n 


n-1 


+ 


c,  h  ) 
1  n 


to  get 


N 

rn-1 

!  k,  T  I I r . I I .  exp 

1  'n  "  1  11 -3  11  3 

|0  kl(col6j+l_0jl+clhj+l). 

N 

<  kx  exp(k1cQV  +  k1c]L) 

thus  proving  stability.  □ 


With  little  additional  effort  it  is  possible  to 
enlarge  the  class  of  problems  for  which  it  is  known 
that  rth  order  convergence  is  possible.  The  next  theorem 
requires  only  that  y ^r+P  ^  (t)  be  of  bounded  variation. 
This  idea  has  been  used  by  Chartres  and  Stepleman  (1972) 
to  analyze  Euler's  method. 

Theorem  5.2: 

Consider  a  multivalue  method  which  satisfies  the 
rth  1-order  condition.  Let  y(t)£Cr+p  [0,1]  be  the 
solution  of  a  Lipschitz  problem  where  y ^r+p  ^  (t)  is 
of  bounded  variation  on  [0,1].  Then  aS  (t ;h)  =  DY (t ;h) 
satisfies 

1.  ||aS(t;h)  -  a(t;h)  |  |H  =  0(hr)  . 

2.  a S ( t ; h )  =  aS(t;h)  +  0 (hr+P_1v [ t-h , t] )  +  0 (hr+P) 

where  v[t-h,t]  is  the  variation  of  y^r+p  1  (t)  on  [t-h,t] , 
which  shows  the  method  is  consistent  of  order  r. 


Proof:  With  one  modification,  the  proof  of  Theorem 

3.3  can  also  be  used  here.  When  Y(t;h)  is  approximated 
by  PY (t-h ; h) ,  the  jth  component  of  the  error  is 


r+p-1 

3 


(y(r+p_1)  (t-h)-y (r+p-1)  (xj) » 
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where  Tj£(t-h,t) 


PY (t-h ; h) 


depends  on  t  and  h.  Therefore 
=  Y  (t ;h)  +  0  (hr+P  1v[t-h,t])  . 


□ 


Theorem  5.3: 

Consider  a  multivalue  method  which  satisfies  the 
1-root  condition  and  the  rth  1-order  condition.  Let 
y  (t)  [0 ,1]  be  the  solution  of  a  Lipschitz  problem 

where  y^r+p_1^  (t)  is  of  bounded  variation  on  [0,1]. 

Then  the  computed  solution  satisfies 


<  kn 
n  -  0 


N 

I 

n=0 


r 

— n 


n 


+  0(hr) 


which  shows  the  method  is  convergent  of  order  r. 


Proof:  Let  ajj  ,a®  , .  .  .  be  the  computed  solution  with 


Eo  =  c'VtChp  -  a(0;h0) 


Hn  =  CniiaS(tn;hn+l>  "  ( W 


where  hN+1  =  hN  and  a  (t;h)  =  DY(t;h)  .  Then 


S  =  C.!,aS(t  ;h  ,,)  . 


— n 


'n+1—  v  n'  n+1 


Using  the  triangle  inequality 


■ 


101 


N 


(5.3)  l 


n=0 


s 

r 

— n 


n 


N 

£  I 

n=0 


Cii+1— S  (tn'hn+X)  S  (tn'hn 


n 


as  (0 ;hQ ) -a (0 ; hQ ) 


0 


N 

I 

n=l 


2r(VV~S.  (W 


n 


Let  us  bound  the  terms  of  the  first  sum  on  the  right-hand 
side : 


(Cn+X^S(tn!hn+l> 


r+p-1  .  .  .  y^  (t  ) 

l  |  (cT^h -  h^)D.. 

j=0 


n+1  n+1  n  i  j  j I 


,Vj>  (t„) 


s  le  o  I-T1  A~lj~iLl  D  ^ _ r-L 

1  1 0n+l  6n'  1  -  A  Dij  j  l 


The  first  min(i,p-l)+r  terms  of  the  sum  are  zero;  hence 


Cn+1— S  ^tn'hn+1)  -  aS(tn;hn)||n  =  0(h^|6n+1  -  0n|) 


From  Theorem  5 . 2 


S-S(tn;hn)  -  ^S(tn;hn 


n  =  0(hnv[tn-l'tn1)  +  0(hn+1) 


Clearly 


(5.4)  l|aS(tn;hn)  -  2.<  VV  |  |n  =  O  (h£) 


A 


'  i  *  3' 
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Put  these  last  three  equations  into  (5.3): 


N  N 

y  I  I  r  I  I  =  y  0  (hr  I  6  , 

^  1  1  — n 11 n  L  n  n  1  n+1 


n=0 


n=0 


0nl)  +  o(hp 


N 


r+1 


+  l  (0(hnv[tn_lftnl)  +  0(hn  )) 
n=l 


(5.5)  =  O (hrV)  +  0  (hr )  +  0 (hrv [0,1] )  +  0(hr) 


From  (5.4),  (5. 5), and  Theorem  5.1,  it  follows  that 


an-a(tn;hn) 


1  1  < 

II  x  I 

a  -a 

1  “f"  1 

I  S  X  I 

a  -a 

!  +  1  1 

1  ln  - 

11  — n  — n  1 

1  n  1 

i  ._n  — n  i 

1  n  11 

a  -a (t  ;h  ) 

— n  —  n  n  '  1  n 


<  k 


N 

n=0 


rj|n  +  k0  l  ||j£||n 

n=0 


+1 |H 

1 1  n  n+1  n+1 


afa (t  ;h  ,  ) 
—  n'  n+1 


-(tn;hn+l)  I  I n+1 


=  k 


N 

o  Zo 


r  II  +0 (hr)  +  0 (h^) 

— n  1  1  n  n 


thus  proving  the  theorem. 


□ 


It  has  been  shown  that  the  1-root  condition  and 
the  rth  1-order  condition  are  sufficient  for  convergence 
of  order  r.  That  they  are  also  necessary  follows  from 
Chapter  2  because  fixed  stepsize  methods  are  a  special 
case  of  variable  stepsize  methods. 


. 
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Section  5.2  Adams-Bashforth-Moulton  Methods 

For  these  methods,  p  =  1,  A  =  A,  1^=1,  and 

are  chosen  so  that  the  extraneous  eigen¬ 
values  of  S  are  all  zero.  Therefore  we  expect  excellent 
stability  properties.  The  relative  stepsizes  0n  are 
required  to  satisfy  A  <  0n  £  1  for  some  constant 
0  <  A  £  1  independent  of  h. 

Tu  (1972)  shows  that  an  Adams  method  is  stable 

if  the  value  of  h  does  not  change  more  often  than 

n 

every  k  steps.  Theorem  5.5  improves  this  result  by 
showing  that  it  is  sufficient  to  require  only  k-1 
steps  between  changes  in  stepsize. 


Lemma  5.1: 

Consider  an  A-B-M  method  with  stability  matrix 

S  and  a  step  selection  technique  with  the  property  that 
n-  j 

TT  SC  .  .  is  uniformly  bounded  for  all  j,  n,  and  N 

.  ,  n+l-i 
i=l 

such  that  0  £  j  £  n  <  N .  Then  the  method  is  convergent 
of  order  r  where  r  =  k  except  for  the  2-value  trape¬ 
zoidal  method  UQ  =  j)  in  which  case  r  =  2. 


Proof:  To  prove  stability,  proceed  as  in  Theorem  5.1 

with  a  vew  modifications: 

Write 

e  =  SC  e  n+£he  +  r 
-in  n-n-1  -  n  n  -n 


■ 
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and 


1  e 

<  l| 

|C  I 

1  e 

1  n 1 

1  n  1  1 

i  i  _n 

Solving  the  difference  equation  gives 


n  m-  j 

e  =  y  n  sc 

— n  L 


j=l ^i=l 


n+l-i 


n  rn-i 


£h  .  e  .  +  l 


:  D 


n  sc 


j  =  0  ^i=l 


n+l-i 


and  by  the  hypothesis,  the  matrix  products  are 
bounded  by  some  constant. 

The  A-B-M  methods  have  the  property  that  if  y(t)eC  [0,1] 
is  the  solution  of  a  Lipschitz  problem  then 


y*  J_  "I 

a(t;h)  =  a(t;h)  +  0(h  ) 

where  r  is  defined  above.  The  remainder  of  the  proof  is 
similar  to  the  proof  of  Theorem  5.3. 


Theorem  5.4; 

If  k  =  1  then  the  A-B-M  method  is  convergent  of 
order  1  for  SiQ  ?  j  and  of  order  2  for  =  j. 


Proof : 


n-j 

n 

i=l 


SC 


n+l-i 


1 

0 


0 


if  n  >  j+1 


□ 


* 
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Theorem  5.5: 

Let  k  >  2,  and  assume  that  at  least  k-1  steps 
are  taken  between  changes  in  stepsize.  Then  the  result¬ 
ing  A-B-M  method  is  convergent  of  order  k. 


Proof:  S  satisfies  its  characteristic  equation 

gk+l  _  gk  _  Thus 

S(Sk6  .)  =  Sk6  . 

-3 


and  so  Sk5_^  is  an  eigenvector  of  S  corresponding  to  1. 
The  eigenvalue  1  has  multiplicity  one  and  S6_Q  =  . 

It  follows  that  Sk6_ .  is  a  scalar  multiple  of  6^  for 
j  =  0,1,..., k.  Hence  only  row  zero  of  S  is  nonzero. 


Thus 


Ao^sk 


r 


which  implies 


and 


k-1  ,  „  „T,  r.  rT _k~-l 

S  (I  -  £  o_±)  =  6_q 6 QwS  A 


S1"'1  =  sk-h 


Since  6_^C^S  =  — ' ^  ' 


v _ i  t  k~  — 1 

Sk  "C.S  =  S^S  a  ciS  . 


(5.6) 


The  theorem  follows  from  Lemma  5.1  if  it  is  shown  that 


(5.7) 


n-j 

n 

i=l 


SC 


n+l-i 


is  uniformly  bounded  for  all  j,  n,  and  N.  There  are  two 
possibilities  to  consider. 

Case  It  n-j  <  2k-2.  Since  A  <  an  £  A-1,  each  of  the 
matrices  in  (5.7)  can  be  bounded,  and  so  the  product  is 
bounded. 


Case  II:  n-j  >  2k-2.  If  Cj  +  2  =  ...  =  Cj+k  =  I  then 

define  1  =  2;  otherwise,  choose  l  so  that  2  <  £<  k  and 

C.  n  ±  I.  In  either  case  it  follows  that 
j  +  £  ' 

Cj+£+l  =  •**  Cj+£+k-2  =  1 

whence  the  product  (5.7)  becomes 


n-j -£-k-2 

n  SC 

i=l 


n+l-i 


qk  — 1„  i/-. 

s  cj+r  cj+i  * 


From  (5.6)  we  get 


sk~lcj+j>sil  lcj+i“  ^o^osk^  lcj+£s  lcj+i 


Notice  that  SC^  =  6_q  ,  and  therefore 
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n-j 

n 

i=l 


SC 


n+l-i 


(60^)Skr1Cj  +  JlS^1Cj  +  1 


The  right-hand  side  is  the  product  of  at  most  2k+3 
matrices  each  of  which  can  be  bounded. 


If  there  are  only  k— 2  steps  between  changes  in 

stepsize,  then  the  method  is  no  longer  stable.  For 

example ,  Tu  (1972)  has  shown  that  for  k  =  3  the  eigen— 
n-j 

values  of  II  SCn+^_^  are 
i=l 

n  a?  (a  .  -  1) 

i,0f0,  n  2 

i=  j+1 


For  the  sequence  of  stepsizes  h ,h/10 ,h ,h/10 , . . .  the 
method  is  unstable  because  the  last  eigenvalue  is 


unbounded. 


Section  5.3  g-conyergence  and  Concluding  Remarks 

Consider  the  use  of  a  q-convergent  method  with 
fixed  stepsize  to  integrate  an  appropriate  Lipschitz 
equation.  The  computed  value  a^^  ^ s  a  startin9  value 
for  an  integration  on  [1/2,1],  but  if  q  >  1  then  i!h 
may  not  converge  to  y^  (1/2)  for  p-q  <  i  <  p-1.  Even 
though  it  seems  that  some  vital  information  may  be  lost, 
the  information  is  merely  transformed,  and  the  method 


recovers  the  initial  values  at  t  =  1/2  by  taking  certain 

linear  combinations  of  the  components  of  a^^*  However, 

if  a  different  q-convergent  formula  is  used  for  the 

interval  [1/2,1],  then  this  information  is  not  recovered 

and  q-convergence  is  not  assured  on  that  interval. 

Both  formula  changing  and  stepsize  changing 

operate  on  the  assumption  that  approximates  a(tn) 

and  not  aS (t  ) .  Therefore  in  the  case  of  an  rth  order 
—  n 

q-convergent  method,  these  two  operations  generally 

introduce  an  error  of  order  hr+^*  ^  in  the  computed  value 

a  ,  and  if  a  is  regarded  as  a  new  starting  value,  then 
— n  — n 

it  is  reasonable  to  expect  nothing  better  than  1— conver¬ 
gence  of  order  r-q+1  on  the  rest  of  the  interval. 

A  close  examination  of  Section  5.1  convinces 
one  that  it  is  difficult  to  prove  variable  stepsize 
q— convergence  unless  one  assumes  the  1— root  condition, 
part  1  of  the  1st  1-order  condition,  and  part  3  of  the 
1st  q-order  condition.  It  is  expected  that  there  are 
few  methods  which  are  variable  stepsize  q— convergent 
and  yet  are  not  1-convergent  for  q-dif ferential  equations 
Furthermore,  it  is  difficult  to  obtain  necessary  and 
sufficient  conditions  for  variable  stepsize  q-convergence 
Therefore  it  is  suggested  that  the  concept  of  variable 
stepsize  q-convergence  should  be  abandoned.  That  leaves 
us  with  fixed  stepsize  q-convergence.  But  fixed  step- 
size  multivalue  methods  are  of  little  practical  value 
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because  the  corresponding  multistep  methods  are  more 
efficient.  Therefore  it  is  not  worthwhile  to  study 
q-convergence  with  q  >  1  for  multivalue  methods. 

The  importance  of  multivalue  methods  is  due  to 
the  convenience  of  the  scaled  derivative  representa¬ 
tion.  It  cannot  be  said  that  multivalue  formulas  are 
a  generalization  of  multistep  formulas  because  it  was 
noted  in  Section  2.1  that  the  numerical  solution  of  a 
convergent  corrector  only  method  is  also  the  solution 
of  a  stable  multistep  method.  Thus  multivalue  formulas 
possess  all  of  the  limitations  of  multistep  formulas. 
For  example,  Gear  (1973)  shows  that  the  order  of  an  A— 
stable  multivalue  method  cannot  exceed  two. 

The  analysis  of  this  thesis  has  been  limited  to 

M 

corrector  only  methods  and  P(EC)  methods.  However, 
the  theory  can  be  extended  to  include  Nordsieck-type 
methods  which  use  off-step  points  or  derivatives  of  f. 
The  formula  for  both  kinds  of  methods  can  be  written 


a 

— n 


=  Sa  , 
— n-1 


+  —r  $ (t  ,a  ,  ;h) 
p!  —  n  — n-1 


If  f  =  f  (t) ,  then  there  is  some  operator  ft,  whose  com¬ 
ponents  are  power  series  in  hD,  such  that 


$  (t  ,a  ,  ,*h)  =  ftf  (t  )  . 

—  n'— n-1  —  n 


The  optimal  interpretation  A*  is  then  given  by 


' 


A* 


Se“hDA*  +  £2 


(hD)p 

Pi 


To  show  that  a  method  is  rth  order  consistent,  it  is 
sufficient  to  show  that 


|  |  (A*-A)y(t)  |  |H  =  0(hr)  , 

hp|  |$(t,Ay(t-h)  ;h)  -  £y(p)  (t)  |  |H  =  0(hr+1)  . 


And  to  show  stability,  it  is  enough  to  show  that 
| | H—1SnH | |  is  bounded  and  that  £  satisfies  a  Lipschitz 
condition  in  its  second  argument. 
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